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I. INTRODUCTION 


I. 1. Background and Objectives 

At least three classes of feed system instabilities have been 
identified for liquid propellant rockets 1, 2, 3, 4 ^ each case, the 

dynamic behavior of the feed line plays a major role in the instability. 
Probably the best known problem^ ^ involves a closed-loop coupling 
of the vehicle longitudinal structural modes, the feed system and the 
engine. A second problem^ showed up on Saturn flight AS-504 when 
the S-II center engine displayed severe pressure and thrust oscilla- 
tions rather late in the stage burn. In this case, the problem was 
localized, involving a closed-loop instability of the engine support 
structure, the feed system, and the engine. A third type of problem 
which can exist^ involves only the feed system, the inducer-pump 
combination and the engine. This instability is traceable directly to 
the inducer dynamics, with the feed system and engine providing an 
impedance loading which governs the resultant oscillation frequency. 

All of the above instabilities pose a threat to the success of a 
mission. For simulation purposes, it is readily recognized that ade- 
quate modeling procedures are needed for each portion of the vehicle 
involved in a given instability. The objective of this study was to de- 
velop an analytical method, and its corresponding computer program, 
to allow the study of disturbances of liquid propellants in typical engine 
feedline systems. This method was to include (1) the effect of steady 
flow, (2) the influence of distributed compliances and axial wall stiff- 
ness, (3) the effects of local compliances, such as vapor bubbles, and 
(4) various factors causing coupled responses between the line and 
structure, i. e. , bends, mounting stiffness, forced changes in line 
length, and bellows and PVC couplings. 

In all respects, this objective has been satisfied. This report 
summarizes the modeling results for each item defined above, and 
presents a generalized computer program which allows the user to 
easily determine dynamic performance of a line containing any number 
of elemental components. Use of this computer program is illustrated 
by way of a detailed set of instructions contained in Appendix H, plus 
a presentation of several example problems. 


Superscript numbers in the text refer to References presented on 
page 52 of this report. 
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I. 2. 


General Feedline Problem 


The problem of modeling a given feedline can best be illus- 
trated by an example. Figure 1 shows a hypothetical case involving 
a propellant tank, a feedline, a pump and an engine. In general, 
each element of the model is elastically supported and can experi- 
ence some motion relative to a frame of reference on the vehicle. 
These accelerations, plus pressure oscillations in the engine com- 
bustion chamber and cavitating inducer instabilities, represent dis- 
turbance inputs to the feed system which can cause flow and pres- 
sure oscillations. These flow disturbances may, in turn, couple 
back with the engine and/or structure to produce a large amplitude, 
closed- loop instability. The modeling problem, with respect to the 
feed system, is to describe analytically the dynamic pressure and 
flow at an arbitrary point in the system, with proper treatment of 
the end conditions, the acceleration inputs and the flexible supports. 

Generally, a feedline may be regarded as a series of elements, 
such as lengths of line, bends, bellows PVC coupling, etc. There- 
fore, it is necessary to describe each element individually, and then 
define means for mathematically combining these elements so as to 
produce the proper overall system mode. In this report, each ele- 
ment is basically defined by a four- terminal representation in the 
Laplace domain. Overall system model formulation is also achieved 
in the Laplace or operator domain. Where problem solution is re- 
quired only in the frequency domain, this Laplace formulation is the 
one to use, and the present computer program is set up on this basis. 
For problem solution in the time domain, the Laplace operator for- 
mulation must be converted to an equivalent time domain formulation 
(rational approximate model). The procedure for accomplishing 
this conversion is described herein. 
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3447 

Figure 1. Hypothetical Feed System And Engine 
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II. GENERAL LINE MODEL WITH MEAN FLOW 


II. 1 . Physical Picture and Requirements 

An illustration of the physical picture of the flow behavior in 
a feedline is given in Figure 2. For the present, only the fluid dy- 
namic effects in a line, assumed cylindrical, straight, and with in- 
finitely stiff walls, will be considered. 

The line shown in Figure 2 is assumed to contain a turbulent- 
ly flowing viscous liquid. The mean vector velocity is represented 
by v 0 , which contains only the axial component, v 0 (r). Superim- 
posed on this mean flow are two dynamic velocity perturbation 
quantities, v and v t . The velocity component v represents the 
dynamic velocity perturbation of interest, and vt represents the 
turbulent velocity fluctuations. Similarly, the total pressure is the 
sum of a mean quantity, p 0 , a dynamic perturbation part, p, and 
a turbulent part, p t . 

By substitution of the assumed summation forms for velocity 
into the Navier- Stokes equations^, and short-time averaging of the 
turbulence quantities, the following form results: 

|^+V 0 |^ = - ^ + V {f v ( v * v) - VX(7Xv)}+ TDT (1) 


In addition, we may complete our description of the dynamic 
fluid behavior with a continuity equation: 



- dp 

Po v - v + v ° = 


o 


( 2 ) 


and a liquid equation of state: 

dp = (3) 

Equation (1) is a first-order form of the Navier-Stokes equations, 
including turbulent dissipation terms (TDT). Similarly, Equation 
(2) is a first-order continuity expression relating the dynamic velo- 
city perturbation to density perturbations (p 0 is the mean density). 

Examination of Equations (1) and (2) shows that the presence 
of a mean turbulent flow has two possible effects. First, the 
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TURBULENT PROFILE 



v T ( r,z,t ) = vector total velocity 

V V v+v t 

v 0 =tlme invariant mean flow 
v= dynamic perturbation of interest 
v t = turbulent fluctuations 

P T = Po + P + Pt 

2747 

Figure 2. General Physical Picture Of Flow Behavior In Feed Line 
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turbulent dissipation terms introduce damping in addition to the con- 
ventional "laminar" damping accounted for by the terms in brackets 
on the right side of Equation (1). Second, the mean flow introduces 
a "convective" effect represented by the terms v 0 Sv/dz and 
v 0 Sp/dz in Equations (1) and (2) respectively. 

The requirements for a dynamic feedline model with a mean 
flow are that (1) the attenuation effects of the turbulent flow be pro- 
perly treated and (2) the possible convective influences be examined 
and included if necessary. The following sections treat these two 
topics. In general, the feedline model presented will be an exten- 
sion of an existing distributed parameter model which exactly ac- 
counts for laminar viscous effects, but negligible mean flow. This 
existing model, which is discussed in References 5 and 6, may be 
given in the following form: 

P 2 (s) = Pi (s) cosh YE - Z o A 1 Qj (s) sinh YL. 

Q 2 (s) = Qi(s) cosh yL - Z 0 1 APj t (s) sinh yL (4) 


The laminar propagation operator takes the familiar form 


Y lam ~ c t 1 


2Ji(gr 0 ) I 2 
5r o J 0 ( §r 0 ) - 


(5) 


where c is the phase velocity in the fluid, J 0 and J x are zero and 
first-order Bessel functions of the first kind and 

= -s/v (6) 

The characteristic impedance is related to the propagation operator 
by 


Z c - Po c o Y/s (7) 

The magnitude of the additional attenuation due to turbulence 
is a function of the Reynold's number and the frequency range of 
interest; this aspect is discussed in Section II. 3. The relationship 
between mean flow velocity and the adiabatic speed of sound dictates 
the phase velocity to be used in Equation (5). Section II. 2 discusses 
the effect of mean flow convection. 
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II. 2 


Possible Convective Influences 


As indicated previously, the terms v 0 dv/dz and v 0 dp/dz 
in Equations (1) and (2), respectively, represent the convective 
effect of the mean flow v 0 . The question to be answered at this 
point is: Are these terms of enough significance to be retained in 
the final model? Ideally we would like to retain the convective terms 
at all times, but, from a practical standpoint, they increase con- 
siderably the complexity of the final solution. 


To judge the importance of the convective terms, it is con- 
venient to consider a nondissipative, one -dimensional form of Equa- 
tions (1) and (2) or 


dy z dy z 1 dp 

dt + V ° dz p 0 dz 


and 


dp d Vz 

aT + p " Tt 


+ V 0 


dp 

dz 


0 


( 8 ) 


(9) 


Transformation of Equations (3), (8), and (9) into the Laplace domain, 
and the subsequent elimination of p and P, yields one equation in 
V z (V z is the transform of v z ), or 


(Co' 


Vo 2 ) 


dfVj 

dz 2 


o dV Z 

2s v ° it 


s 2 V z 


= 0 


( 10 ) 


where 


Co = (H/Po ) 2 


Laplace variable 
isentropic speed of sound 


V 0 


= axial mean flow velocity 


Assuming a solution for Equation (10) of the form 
Yz 

Vz = Ce Y 


(id 


we find the characteristic equation to be 

(Co 2 - V 0 2 ) Y 2 - 2s v 0 Y - s 2 = 0 (12) 

whose roots or values for y are 
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( 13 ) 


iV 



where c* 2 = cf - v 2 . The plus sign applies to waves traveling in 
the retrograde direction, and the negative sign to waves traveling in 
the mean flow direction. 


From these results, it is found that for propagation in the +z 
direction, with the mean flow, the phase velocity is greater than c 0 . 
To demonstrate this result, recall that the propagation operator is 
functionally equivalent to the Laplace variable divided by a velocity 
quantity. Substituting the equivalent expression for c* into Equa- 
tion (13), the following result is obtained: 


Y 


s (v 0 ± Cp) 

(v 0 + Co) (v 0 - Co) 


(14) 


For propagation in the mean flow direction. 


Y 


+z 


s 

Co + V 0 


(15) 


Therefore, the phase velocity for downstream propagation is greater 
than the adiabatic speed of sound. Similar reasoning verifies that 
the phase velocity for retrograde propagation is less than the adia- 
batic speed of sound. 

It is important to point out that for normal conditions of 
liquid propellants flowing in a feedline at velocities up to, say, 50 
fps, the phase velocity with mean flow is so minutely changed that 
this convective effect may be neglected. The reason for this is that 
the mean velocity is much smaller than the liquid sonic velocity, c 0 . 
The only condition under which this might change would be where 
sufficient entrained gas is present in the fluid to cause the new ef- 
fective c 0 value to be, say, an order of magnitude lower than for 
the pure liquid. Then the convective effect should be included. For 
example (see Section III. 1), if the level of entrained gas were to 
approach 0. 1 percent by weight of the total liquid mass, then the 
equivalent speed of sound of the mixture would be approximately 
100 fps. In this case, the speed of sound is only twice as large as 
the maximum mean flow velocity, and the phase velocity would be 
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100 ±50 fps, depending on the propagation direction. This example, 
however, is unrealistic because the gas concentration rarely ap- 
proaches a level where the mean flow velocity is but a small frac- 
tion of the mixture speed of sound. Therefore, for all practical 
purposes, the phase velocity can be taken to be the adiabatic speed 
of sound. 


II. 3 Effect of Turbulence on Damping and Phase Velocity 


It is known that the four terminal Laplace domain represen- 
tation for pressure and flow in a line is extremely accurate when the 
flow is laminar and the propagation operator is of the following form: 


r 


L 


/1L\ 1 


Uo ) 1 

§r 0 Jo(5ro) J 


- Yiam k 


(16) 


In reality, a turbulent flow condition exists for many propellant feed- 
line operating points, as indicated by Reynold's numbers of the order 
of 10 5 to 5 X 10 7 . Therefore, the effect of turbulence must be re- 
flected in the propagation operator. As shown in Figures 3 and 4 
(from Reference 7), the predominant effect of turbulence is to in- 
crease the spatial attenuation at low frequencies. At high frequen- 
cies, the laminar and turbulent attenuations coincide. There is some 
tendency for turbulence to reduce the phase velocity at very low fre- 
quencies and high Reynold's numbers. However, after a thorough 
investigation into existing feedline geometries, propellant properties, 
mean flow velocities and frequency range for structural-hydraulic 
coupling, it was concluded that the effect of turbulence on the phase 
velocity could be neglected. 


For the exact distributed parameter line model, the effect of 
turbulent attenuation was modeled by adding to the existing laminar 
propagation operator, T^, an additional attenuation contribution, 
r 0 , due to turbulence. 

+ r 0 = (r Lr + r 0 ) + ir Li (17) 

The form of turbulence induced increment in attenuation was derived 
from the form of the attenuation factor of a constant I-R-C, lumped 
parameter, representation of a line. That is. 


r 0 = ^Cs \/ IR/s + 1 L 


(18) 


9 




Figure 3. Predictions of Attenuation Factor for Turbulent Flow ( Reference 7 ) 
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Figure 4. Predicted Phase Velocity with Turbulent 

Flow ( Reference 7 ) 


n 



By analogy. 



( 19 ) 


where 


Rt 


2VK N R n 


( 20 ) 


The constants, K and n, were determined by surface fitting the 
turbulent attenuation increment in Figure 3 as a function of Reynold's 
number and frequency. The derived values of these constants are: 


K = 0. 0055 
n =0.85 
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III. DISTRIBUTED COMPLIANCES 


/ 

III. 1 Distributed Bubbles 


In most practical applications, the liquid rocket propellant in 
the feedline will contain quantities of either dissolved ullage gas or 
propellant vapor, or possibly both. It is assumed that the gases or 
vapor are in the form of small bubbles, homogeneously distributed 
throughout the liquid. The net effect of distributed bubbles in the 
liquid is to modify (reduce) the effective fluid sonic velocity and 
possibly to introduce added damping because of irreversible com- 
pression and expansion of the bubbles caused by a periodic distur- 
bance. 


Two types of entrained gases have been considered in the 
model for the acoustic velocity in a two-phase flow: 

(1) The flow of a single- component, two-phase mixture, 
such as LOX and oxygen vapor, and 

(2) The flow of a two- component, two-phase mixture, 
such as LOX with entrained helium ullage gas. 

For a single -component, two-phase mixture, the actual flow 
process is theoretically bounded by the equilibrium process and the 
constant quality process. In the case of an equilibrium process, the 
vaporization and condensation rates are large enough to ensure that 
the vapor temperature and liquid temperature are always equal; 
both phases are saturated, and the quality (vapor fraction by weight) 
varies only with pressure. Based on the assumption of thermodynamic 
equilibrium, Gouse and Brown^ have shown that the speed of sound in 
a single -component, two-phase flow can be expressed as 

C ° = (^77 + X ) (S e- Sf) MH-) v ] (21) 

where x is the mixture quality and (S g - Sf) is the change in entropy 
between the saturated vapor and the saturated liquid state. Evalua- 
tion of Equation (21) requires the extensive use of a temperature- 
entropy chart or a tabular representation of the thermodynamic pro- 
perties of the fluid being considered. 

On the other hand, the condensation and vaporization rates in 
a constant quality process are assumed to be so low that no significant 
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phase change occurs. In a two- component mixture, a constant qual- 
ity exists at all times if there is no liquid vapor present. Rearrang- 
ing the derivation by Hsieh and Plesset^, the speed of sound in a 
two- component mixture is 


1 


cf 


P 


{( 


£i 

Ps + 0P£ 


) 


n 


1 / \ 1 11 

cf + VPg + tfPj i ) P 6 c 8 2 J 


where 


P 

P 8 

Pj& 

and 0 


mixture density 
gas density 
liquid density 

speed of sound in the liquid 
speed of sound in the gas 
mass gas/mass liquid. 


( 22 ) 


As the mass ratio, 0, approaches zero, the speed of sound ap- 
proaches that for the pure liquid. Plesset and Hsieh^O have also 
shown that, for very small gas bubbles dissolved in a liquid, the 
adiabatic speed of sound in the gas should be replaced by the iso- 
thermal speed of sound. The reason is that the gas bubbles have a 
non-uniform temperature distribution in their interior due to the 
finite value of the heat conduction rate. Figure 5 presents the iso- 
thermal versus the adiabatic model for dissolved helium ullage gas 
in LOX. For comparison, the equilibrium model for LOX with 0 2 
vapor is also shown. Reliable experimental data are not available 
to determine which flow process actually occurs for the flow of a 
single- component, two-phase mixture. As can be seen by Figure 5, 
there is a great difference in the acoustic velocity between the equi- 
librium and the constant quality flow processes. These differences 
could lead to a large error in the predicted feedline system re- 
sponse, and more work should be devoted to a study of the actual 
flow process. In the present study, the constant quality model is 
applied to both single- component and two- component flows. A more 
detailed discussion of the acoustic velocity in two-phase flow is pre- 
sented in Appendix B. 


The wave propagation characteristics of the mixture depend 
not only on the gas-to-liquid mass ratio but also on the bubble size 
and disturbance frequency. Bubbles of a given size resonate when 
excited at their natural frequency; this natural or resonant frequency 
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Figure 5. Variation Of Speed Of Sound With Mass Ratio 
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increases in inverse proportion to the bubble radius**. An investi- 
gation has been made to determine the variation in effective sonic 
velocity in a liquid with uniformly distributed gas bubbles as the 
excitation frequency approaches the bubble resonant frequency. 


At frequencies approaching resonance, Equation (22) is no 
longer useful, even as an approximation. Spitzer*^ derived an ex- 
pression for the phase velocity at all frequencies by assuming that 
the gas concentration is small. If surface tension is neglected and 
if the bubble diameter is assumed to be larger than 0.001 ft, the 
bubble resonant frequency can be approximated by*^» *^ 


1 / 3T|P 0 4|a / 

° “ 2ttR 0 y p A ' p/R® ) 


(23) 


where r), the polytropic exponent, is dependent primarily on the 
bubble size and gas thermal conductivity. For very small bubbles, 
such as those assumed to be distributed in the propellant feedline, 
the polytropic exponent is equal to unity. 

The theoretical effect of frequency and bubble size on the 
phase velocity in a mixture of air bubbles in water has been pro- 
grammed and is presented in Figure 6. The mass ratio, 0, was 
assumed to be 3.0 X 10 . Note that the speed of sound is only af- 

fected at or near the bubble resonant frequency, which increases 
with decreasing bubble size. For the homogeneous distribution of 
very small bubbles entrained in the liquid propellant, the net effect 
over the frequency range of interest is only to lower the phase velo- 
city of the pure liquid, to the value predicted by Equation (22). The 
effect of larger bubbles, which are assumed to exist only locally 
(cavitation zones), is discussed in Section IV. 


III. 2 Distributed Wall Compliance 

In general, the flexibility of the feedline wall will produce two 
possible effects which may be of concern. First, the wall compli- 
ance will reduce the phase velocity and increase the spatial attenua- 
tion for the longitudinal fluid wave propagation mode (called the 
zeroth mode in Reference 5), and second, the axial wall stiffness 
effect, will permit real wave propagation of the higher order modes 
at low frequency by virtue of a coupling between the fluid and axial 
wall modes. For typical feedline problems, we have concluded that 
the amount of energy which is fed into the higher order modes at the 
line terminations is so small that it can be neglected. Therefore, 
for all practical purposes, the effect of the axial wall stiffness on 
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Mass Ratio, 4 > = 3.0 x 10~ 5 or 
Volume Ratio, a = .024 
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Figure 6. Effect Of Frequency And Bubble Size On Phase Velocity Of Air And 

Water Mixture 




the wave propagation can be neglected. This topic is thoroughly- 
discussed in Section III. 3. 

The influence of the wall compliance (and mass) on the pro- 
pagation of the zeroth mode cannot be neglected. Gerlach^ investi- 
gated the case of a viscous fluid flowing through a line having elastic 
flexible walls, i. e. , a line whose walls have finite radial compli- 
ance and impedance but do not propagate disturbances axially. By 
equating the radial impedance of the wall to the radial impedance of 
the fluid, the latter being subject to the no- slip condition on axial 
fluid velocity at the wall, the following characteristic equation for 
the eigenvalue problem was obtained: 

(V 3 - P 2 ) Pocf 

[p Ji(Pr 0 ) Y» Ji(kr,) l = p t hs s + hE t /r? (24) 

L J 0 (Pr 0 ) " k J 0 (kr 0 ) J 

The parameters, |3 and k, are separation constants of the solution 
of the Navier-Stokes equations and are related by 

Y 2 = k 2 + s/v = P 2 + s 2 /c 0 (25) 


where 


h 

r 0 

E t 

Pt 

s 


C o 
V 


Po 

j 0 » Ji 


tube wall thickness 
tube radius 

Young's modulus for tube material 
density of tube wall 
Laplace variable 

isentropic speed of sound in the fluid 
fluid kinematic viscosity 
fluid density 

zero and first-order Bessel functions of the 
first kind. 


Equations (24) and (25) were solved numerically on the computer. 
Figures 7 and 8 illustrate the effect of an "elastic, flexible wall" 
(a wall with radial compliance and mass) on the spatial attenuation 
and phase velocity for the zeroth mode. The area of interest in 
these figures corresponds to dimensionless frequency numbers, 
DJr 0 /c 0 , less than 1.0, since it is known that coupled structural- 
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Figure 7. Zeroth Mode Spatial Attenuation Versus Frequency 
Number For A Rigid And An Elastic Flexible Wall 




Figure 8. Zeroth Mode Dimensionless Phase Velocity, c /c 0 , Versus 
Frequency Number For A Rigid And An Elastic Flexible Wall 
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feedsystem instabilities occur in this frequency range. In this range, 
the wall behaves in a spring-like manner and inertia effects are neg- 
ligible; radial attenuation for the elastic wall is virtually indistin- 
quishable from the attenuation of the rigid wall. It has been estab- 
lished that for any practical feedline problem, this is the case; 
hence, only the wall radial compliance need be considered. Wall 
mass effects are negligible. Therefore, the classic Korteweg cor- 
rection to the phase velocity is valid. That is, 

c = F#F 

It must be remembered, however, that this correction is the result 
of applying a boundary condition to the fluid dynamic formulation, 

III, 3 Axial Wall Stiffness 


In the preceding section, it was stated that the effect of axial 
wall stiffness on higher mode wave propagation could be neglected. 
That statement can be justified by proving the following points: 

(1) Due to observable flow processes in feedlines, there is 

a minimal amount of energy available in the higher modes 
of the fluid for coupling with the wall, 

(2) The energy that is available in each higher fluid mode 
has an attenuation rate, below the cut-off frequency, of 
approximately two orders of magnitude larger than that 
for the zeroth mode. Therefore, the energy level and 
attenuation rate do not permit sustained coupling of the 
fluid and the wall in the higher modes. 

The physical effect of higher modes must be considered, because the 
proposed distributed parameter line model is based on only the 
zeroth mode of propagation, and the existence of substantially 
higher modes would invalidate its functional form. 

Recall that the magnitude of Reynold's numbers in a typical 
feedline indicates a turbulent flow condition. Hence, the velocity 
profile should be relatively flat. Gerlach^ has calculated the pertur- 
bation velocity profiles of the first four modes for the laminar, 
viscous pipe-flow problem. Characteristic results are shown in 
Figures 9 and 10 where F Zn (r) is the velocity distribution function. 
Observe closely that, for the wide range of frequencies and damping 
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Figure 9. Axial Velocity Profile Function, F zn (r), for Four 

Modes (ajr 0 /c 0 = 1.0, Wr 0 c 0 = 0.01 ) 
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Figure 10. Axial Velocity Profile Function , F zn ( r ) , for Four 

Modes (wr 0 / c 0 = 0.2, i//r 0 c 0 * 0.001 ) 
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numbers considered, the zeroth mode velocity distribution is near- 
ly constant across the cross-section. The obvious conclusion is 
that the zeroth mode perturbation profile approximates a turbulent 
profile, and, hence, a minimal contribution from the higher modes 
is required to fill out the turbulent profile. Therefore, it follows 
that the energy level in the higher modes is negligible. 

Based on the foregoing discussion, there is admittedly a 
small contribution of the higher modes to the velocity profile. How- 
ever, Gerlach^ has also shown that the spatial attenuation (the real 
part of the propagation operator) for these higher modes is substan- 
tially larger than for the zeroth mode. The net effect is that the 
higher modes are rapidly damped out over a broad frequency range, 
and therefore have neither the time nor the energy levels required 
to excite (couple) corresponding modes of transmission in the con- 
duit wall. Figure 11 illustrates the elevated damping level of higher 
modes, while the zeroth mode attenuation is roughly insensitive to 
the frequency range. 

The above discussion substantiates the validity of the zeroth 
mode transfer equations for the line. In addition, the foregoing con- 
clusions which were derived on a qualitative basis, have also been 
verified quantitatively by Lin and Morgan^. 
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IV. LOCAL COMPLIANCES - LARGE BUBBLES AND 
CAVITATION ZONES 


At elbows or bends in the feedline, or at the pump inlet, local 
pressure drops are likely to create regions of cavitation. In most 
cases, these regions are composed of localized bubbles of vapor. 
These larger cavitation bubbles in the line are modeled by consider- 
ing the bubbles to be local compliances. As the excitation frequency 
approaches the bubble resonant frequency, which can be quite low 
for large bubbles, the volume pulsations may be described by the 
following linear, second- order differential equation: 


PX 

4rrR 0 


v + bv + 



-p(t) 


(27) 


or 


msv + bv + kv = -p(t) 


(28) 


where 

= liquid density 
P 0 = steady state line pressure 
V 0 = steady state bubble volume 

R 0 = steady state bubble radius 

b = coefficient accounting for viscous, thermal and 
radiation damping. 

The continuity equation relating the flow upstream and down- 
stream of the bubble is 


q 2 - qj. = dv/dt (29) 

where qx and qg represent the upstream and downstream volu- 
metric flow rate, and dv/dt is the change in bubble volume with 
time. Equations (28) and (29) can be combined conveniently in the 
Laplace transform domain to yield an expression for the change in 
flow rate due to the presence of the bubble. 


Qa(s) - Qi (s) 


- s P(s) 

(m 2 s 2 + bs + k) 


(30) 
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where the term l/(m 2 s 2 + bs + k) can be visualized as a "complex" 
compliance. The damping term, b, has been estimated by Devin^, 
and is represented by a combination of thermal damping, sound 
radiation damping, and viscous damping. 

For zero damping of the oscillations, 


Qi (s) - Q 2 (s) 


sP(s) 


(31) 


At frequencies considerably below the bubble natural frequency, uu 0 , 
Qi(s) - Q.(.) (32) 

where uu 0 is defined as 

it) o = V k/m 2 ’ (33) 

For these low frequency cases, the compliance becomes simply 


C 



(34) 


The correct value of the polytropic exponent r\ in the stiff- 
ness term k = r\ P 0 /V 0 depends on the bubble size and frequency, 
as well as the thermal properties of the gas and liquid. This varia- 
tion with bubble size and frequency of excitation is shown in Figure 
12. For a given bubble undergoing successive compressions, the 
gas near the center behaves adiabatically, while the gas near the 
liquid-gas interface undergoes no change in temperature since the 
liquid behaves as a large heat sink. The polytropic exponent, T), 
is calculated by the method described in Appendix C. 
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Figure 12. Polytropic Exponent for 0 2 Bubble in LOX vs Dimensionless Frequency and Bubble Size 



V. COUPLED RESPONSES 


V. 1. Accelerated Line 


One possible mode of structural-hydraulic coupling is shown 
in Figure 13. An externally imposed, rigid-body acceleration acts 
on a length of feedline that is transporting liquid propellant. Be- 
cause of viscous shearing forces at the propellant- conduit interface, 
the imposed line motion will couple with the fluid motion to produce 
modifications to the pressure-flow equations for the stationary line 
(Section II. 1). The functional form of these modifications can be 
easily derived. 


The fluid may be characterized by the familiar axisymmetric 
linearized, perturbation equations of motion, continuity and state. 


du 

at 


1 9p rd 2 u 

p 0 Sz + V Ur 2 


+ 


1 

r 



(35) 


ap 

at 


+ 



o 


(36) 


ap _ ap 

Po " K ’ 


H — P o C o 


(37) 


Implicit in this formulation is the assumption that the radial pertur- 
bation velocity is negligible compared to the axial perturbation velo- 
city, u. Note that the pressure has not been assumed to be inde- 
pendent of the radial coordinate. The body force F z , can be re- 
placed by the equivalent D'Alembert force, -p 0 a z , due to the im- 
posed acceleration. Utilizing this substitution, combining the con- 
tinuity and state equations and applying the Laplace transformation 
with zero initial conditions, yields the following pair of coupled 
differential equations: 


sU 



ap 

dz 


+ V 


rafu 

Ldr 2 


+ 


1 _ 

r 



(38) 


PoCq 2 au 

s dz 


(39) 


where P and U represent the transformed fluid pressure and velo 
city. These equations can be combined into a single equation in the 
dependent variable U. 
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sU = 


5 2 U 
dz 2 


ra?u 1 5ui A 
v L Sr 2 + r dr J " Az 


(40) 


Making use of the fact that A* is a function only of the Laplace 
variable, s, and redefining the dependent variable to be 

V (r, z, s) = U + A z /s (41) 

results in the following easily solved, ordinary differential equation: 

v r 5 g v i av _ s ~] c„ a 

s Ldr 3 ^ r dr v J s 2 


2 dfv 

dz 2 


(42) 


This equation is amenable to solution by the standard separation of 
variables technique which yields 


V = (A cosh Xz + B sinh Xz) J 0 (ar) 


(43) 


where 



X 2 = separation constant 
A, B = integration constants. 

By definition, X is the propagation operator, Y, which governs the 
spatial attenuation of axially propagated waves in the fluid. This 
operator may be stated a priori: 


Y 


s 

Co 


( l 


1 

2Ji(gr 0 ) 
?r o J 0 (5r o ) J 


(44) 


where 


2 

5 = -s/v , r 0 = conduit radius 

In addition, the characteristic impedance of the line, Z c , is defined 
to be 


Z c = 


Po c o 


S 


Y 


(45) 
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Utilizing these definitions, the fluid pressure and velocity become 

U (r, z, s) = [A cosh Yz + B sinh Yz] J 0 (ocr) - A z / s (46) 

P (r, z, s) = - Z c [A sinh Yz + B cosh Yz] J 0 (ar) (47) 

The dependence on the radial coordinate may be eliminated by 
averaging these expressions across the cross section. 

U = [A cosh Yz + B sinh Yz] — — (ar 0 ) - A z /s (48) 

cxr 0 

P= -Z 0 [A sinh Yz + B cosh Yz] — — Jx (ar 0 ) (49) 

ar 0 

Constants, A and B, can be evaluated at the origin of the line (z=0) 
by applying the boundary conditions. 

U(o, s) = Ui ; P(o, s) = Pi (50) 

which results in the following expressions: 

P 

U = Ui cosh Yz - sinh Yz - A z /s (1-cosh Yz) (51) 

Z c 

P = Px cosh Yz - Z c Ux sinh Yz ^ sinh Yz (52) 

s 

The desired four-terminal representation for the mean exit pressure 
and flow velocity are obtained by evaluating the above expressions at 
z = L. 

U 2 = U% cosh yL - Pi/Z c sinh yL - V z (1-cosh yL) (53) 

P 2 = P x cosh yL - Z 0 U% sinh yL - V z sinh yL (54) 

where the acceleration has been replaced by its Laplace equivalent 
s V z . These equations can be referred to local volumetric flow rate 
by introducing the line cross-sectional area, A. 

V. 2 Relative Motion of Bellows and PVC Joints 

Bellows 

Bellows are commonly used elements in propulsion feed sys- 
tems, and in their most important application, they serve as a "fix 11 
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for the POGO instability. Bellows vary in their structural com- 
plexity from the simple array of axisymmetric convolutes shown in 
Figure 14a to the configuration in Figure 14b which incorporates a 
gas trap liner. The following discussion is directed toward the 
latter configuration which includes the simple bellows as a special 
or degenerate case. 

When the ends of a bellows execute relative motion, the change 
in volume flow rate between the inlet and exit can be considered to be 
the sum of two separate effects: 

(1) A flow rate change due to the compliance of the gas 
trapped under the liner, and 

(2) An apparent local volume production. 

It should be noted that the first item exists even in the ab- 
sence of relative motion. 

First, the pressure drop across a bellows in the flow direc- 
tion is expressed in the Laplace domain by 

Ps(s) = F(PV z 2 ,G) Pi (s) , (55) 

where F(pV z 2 , G) is a loss factor, dependent on the dynamic pres- 
sure and bellows geometry. Now, the change in flow rate due to 
trapped gas can be written in terms of the bellows inlet pressure. 

q 3 ' - qi' = - c (56) 

or in the transform domain 

Qz(s) - 0/(3)= -Cs P x (s) (57) 

where C is the lumped compliance of the gas and is defined by the 
ratio of the gas volume to the gas bulk modulus of elasticity. 

The apparent volume production can be derived from the con- 
tinuity equation 

q 3 " - qi* = dv/dt (58) 

or in the Laplace domain 

Qa"(s) - Qj^s) = -sVol(s) (59) 
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( )' - quantities associated with compliance 
of trapped gas 


( )" - quantities associated with volume 
production 

b. Bellows With Gas Trap Liner 

3442 


Figure 14. Nomenclature For Relative Motion Of Bellows 
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The volumetric production is related to the axial displacement of 
the bellows ends by 


Vol(s) = K b [X x ( s ) - X s (s)] (60) 

where K b is a "volume change" constant, also dependent on the 
bellows geometry. The resulting flow equation is 

Qa'(s) - Qi"(s) = sK b [X x (s) - X 2 (s)] (6l) 

or in terms of relative velocities 

Q/(s) - Q/(s) = Kb [Vi (s) - Vg (s)] . (62) 

The combined effect, expressed as the sum of Equations (55) and (62) 
produces the following result for the change in flow rate: 

Qs(s) - Qx(s) - CsP 1 (s) + K b [V x (s) - V a (s)] (63) 


Equations (55) and (63) describe completely the four terminal pres- 
sure-flQw relationship for a bellows, 

PVC Joints 


In many applications, the PVC (pressure- volume compensa- 
tor) joint is often used in place of an elementary bellows where the 
local volumetric change effect described above is judged undesirable. 
A typical PVC joint, shown in Figure 15, generally consists of a 
bellows to allow line flexibility (either axially, or to allow angula- 
tion, or both) and a compensation chamber to "absorb" the excess 
volume of fluid created in the line when the bellows is displaced 
axially. In theory, the PVC joint is designed to permit axial length 
changes without introducing a corresponding fluid volume disturbance 
in the line. In practice, however, the PVC joint is not a perfect 
compensation device. 

The 'Volume flow production" of the bellows Q 4 , may be 
written as 

Q*(s) = Q a (s) - Q 3 (s) = K b [V x (s) - V 2 (s)] (64) 

where V A (s), i = 1 or 2, is the transformed velocity associated with 
coordinate X* . Similarly, the ideal volume compensation of the PVC 
may be written as 
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Figure 15. Illustration Of Typical PVC Joint 
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Qa(s) = Qi(s) - Q 3 (s) = K 0 [ Vj, ( s ) - V 2 (s)]. 


(65) 


Continuity requires that 

Q 5 (s) = Qi(s) + Q 4 (s) - Q 2 (s). (66) 

Combining Equations (64), (65), and (66), the resultant flow be- 
comes 


Qs(s) = Qi (s) + (K b -K e ) [Vj, (s) - V 2 (s)] . (67) 

It is normally intended that K b = K 0 so that perfect compen- 
sation is obtained. However, there is a dynamic effect associated 
with the action of the compensation unit, so a more general form of 
Equation (65) is 

Q a(s) = [Vi (s) - V 3 (s)] (68) 

where G(s) might be a simple lag, e. g. , (1 + ts), caused by a com- 
bination of the PVC structure elasticity and fluid compressibility. 
Equation (67) was selected as the PVC flow rate expression in the 
feedline computer code. The pressure drop across the joint was 
assumed to have the same functional form as Equation (55). 

V. 3. Forced Changes in Line Length 

Propellant feedlines are normally supported by the vehicle 
structure at discrete points on the line as opposed to a continuous 
support. A forced change in the length of the line between consecu- 
tive supports occurs when vehicle structural inputs to these supports 
differ in phase and/or magnitude. In an idealized fashion, it can be 
assumed that the length change is a result of discrete velocity inputs 
at the extremeties of the line, as shown in Figure 16. 

The effect that length changes produce on the pressure-flow 
relationships for the liquid in the line can be analyzed from two dif- 
ferent viewpoints. The first approach requires the solution of the 
fluid dynamic equations of motion subject to an inhomogeneous bound- 
ary condition on the spatial variation of axial wall velocity as dic- 
tated by the structural inputs, V x and V s . An alternate approach 
assumes that the dynamic length change can be modeled as the linear 
superposition of (1) the axial vibration of a rigid line and (2) a 
volume production region that reflects the relative motion between 
the ends of the line as indicated in Figure 17. The results of each 
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Figure 16 . Geometry For Forced Changes In Line Length 


analysis technique are summarized below, while specific develop- 
ment details are presented in Appendix D. 

The first method of analysis begins with the Laplace trans- 
form version of the linearized (1) radial and axial equations of 
motion for the fluid, (2) continuity equation, and (3) the equation 
of state. These equations are reduced to one -dimensional form by 
averaging across the feedline cross section. The solution of the 
resulting differential equations for the transformed axial velocity is 
matched to the transformed, inhomogeneous, compatibility condi- 
tion on wall velocity. This boundary condition is obtained by solv- 
ing the undamped wave equation for the wall. The resulting trans- 
formed pressure-flow equations are: 


Q 2 -Z c A sinh YL 


cosh YL -Z c /A sinhYL Pi - Y 2 ^V x (s) oq 

A sinhYL coshYL Qi - Y 2 j Ad- 
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Figure 17. Illustration Of Approximate Modeling Procedure For 

Forced Changes In Line Length 
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In these expressions c K is the longitudinal wave speed in the wall 
and G(s) is the transfer function relating the structural velocities, 
Vi and V s . All other terms have been defined in previous sections. 
It is evident that the effect of forced line length changes is to modify 
the pressure-flow relationship for a line with no external excitation. 
It can be shown that when V, = V 2 the above expressions degenerate 
to 


X‘ 


coshyL -Z c A l sinhYL 


X 


Z c sinh yL 


= 

-l 



+ Vf 


_Q ? 


-AZ 0 sinhYL coshYL 


_Qi 


A(l-cosh yh 
— ^ 


where 


_ 2 1,(50 

§r 0 J o( 0 ) 

This latter result was predicted by Gerlach-^ for rigid body, axial 
vibration of a line. 


In the second modeling approach, the pressure-flow equations 
for the rigid body motion of a line are modified with a volume produc- 
tion element that reflects the relative motion of the ends of the line. 
For a rigid line oscillating with velocity, V, , 


' p * 


cosh YL 

-A 1 Z c sinh YL 


X 


Z c sinh yL 

Q s 'i 


-Z 0 X A sinh YL 

cosh YL 


Qi 

- v i 

A(l-cosh yL 


The rate of volume generation due to relative motion is 


Q/ = A(V 2 - V,) (72) 


Total volumetric flow rate is then the sum of Q 2 ' and Q 2 ". Com- 
bining Equations (71) and (72) yields 


Pa' 


cosh YL 

'A 1 Z Q sinhYL 


X 


Z c sinh yL 

QaJ 

- 





-V, 



- Z c A sinh yL 

cosh Y L 


_ Ql 


A(2-cosh yL-G(s) 


As before, G(s) describes the relationship between the externally 
imposed velocities, V, and V 2 . 


Equations (69) and (73) cannot be compared directly since, in 
the first case, the excitation is coupled to the fluid by a boundary 
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condition and, in the latter case, the coupling takes the form of a 
body force applied to the fluid. 

We believe that the result obtained by treating the problem 
as an inhomogeneous boundary valued problem is the most accurate, 
and it is this expression. Equation (69), that has been incorporated 
into the computer code. 

V. 4. Mounting Stiffness 

Propellant feedlines are anchored to the primary vehicle 
structure by mounting brackets that exhibit varying degrees of 
elasticity and damping. These brackets act as filters which modify 
the structural excitation that is transmitted to the feedline. There 
are two techniques that can be utilized in modeling the effect of 
mounting stiffness: (1) discrete parameter -acceleration, and (2) 

impedance. Both of these methods are summarized in this section; 
detailed developments are contained in Appendix E. The transfer 
equations for the discrete parameter -acceleration technique have 
been programmed as Subroutine Eight in the computer program while 
the impedance approach is programmed in Subroutine Ten. 

Discrete Parameter- Acceleration Technique 

Consider the idealized configuration shown in Figure 18. The 
entire line, which has two bends, is assumed to be infinitely rigid and 
elastically restrained. A single linear spring in parallel with a vis- 
cous damper represents the combined effect of all the discrete 
mounting stiffnesses whose line of action coincides with that of the 
equivalent mounting stiffness. In Section V. 1, an expression was 
derived for the pressure and flow rate in a vertical line that is ex- 
cited by a velocity, V^, applied directly to the line. That expres- 
sion forms the starting point for the analysis of the present configura- 
tion. The analytic steps, which will now be summarized qualitatively, 
are presented quantitatively in Appendix E. Initially, the velocity, 
Vj^(s), is replaced by its Laplace equivalent, a^(s)/s. The problem 
is to express the resultant line acceleration, a^(s), in terms of the 
applied structural acceleration, ai(s), the viscoelastic support 
parameters and the fluid mass contained in the horizontal limbs of 
the feedline. By neglecting the elbow at the exit to the fuel tank, 
this aspect of the problem is easily overcome by viewing the entire 
system as a damped, harmonic oscillator whose mass is equivalent 
to the combined liquid mass in the horizontal limbs. The support is 
excited by ai(s), and the inertial effects in the vertical line seg- 
ment are taken into account by applying the dynamic fluid pressure 
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Figure 18. Model For A Line With Mounting Stiffness 



forces on the projected areas of the elbows to the system mass, 
i. e. , at Stations 2 and 3. Solution in the Laplace domain yields 


a^(s)-a 1 (s)[ JvIs2+gb+k ] + A ( p 3- p s) [ Mg 2 +gb+k ] (74) 

where M is the fluid mass of the oscillator. Equation (74) is then 
substituted into the accelerated line equations, and after simplifica- 
tion, the following result is obtained: 
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The L^ terms are the elements of the original matrix for the 
accelerated line. The remaining terms, i. e. , a/, a* g 7 , and f3 y/ , con- 

tain the system mass, support spring constant and damping coeffi- 
cient. As the support becomes progressively stiffer. Equation (75) 
degenerates into the expression for a simple accelerated line 
(Section V. 1). 

Impedance Technique 

The impedance technique will now be applied to the same line 
configuration that was used in the previous case with the exception 
that the mounting stiffness which consisted of a grounded accelera- 
tion, at, and a viscoelastic support, is replaced by a driving point 
impedance, Z 6 . As before, the starting point for the analysis is 
the four -terminal representation of line that is subject to velocity 
excitation, (see Section V. 1). It is required that the line velo- 

city be expressed in terms of the imposed forces and ultimately in 
terms of the driving point impedance. The equation of motion for 
the three- segment line in the Laplace domain is 

F - A(P 3 -P 2 ) - Ms V i (76) 

The quantity, F, represents the lumped external force applied to 
the line through the impedance element, and it essentially replaces a 
combined effect of the spring-dashpot combination in the previous 
development. Using the definition of driving point impedance, 

F/V a, in conjunction with the above equation of motion, the line 
velocity can be expressed by 
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( 77 ) 


V* 


A(P 3 -P 2 ) 
Z s + Ms 


where Z B is the structural impedance. All other quantities have 
the same definition as before. The resulting transfer matrix for 
line response is obtained by substituting this expression for the line 
velocity into the accelerated line equations in Section V. 1. The re- 
sult is 
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where 


( 78 ) 


AZe 2 sinh Tg 
a = Z B + Ms 

_ A 8 (l-coshr s ) 

P " Z B + Ms 

Lii^ = elements of original accelerated line matrix. 


In the case of a rigid support, the structural impedance ap- 
proaches infinity, and the transfer equation degenerates into the 
familiar four-terminal representation of a simple, stationary line. 
It is interesting to note that even though the line is being excited 
externally, the form of the response equations do not include the 
trailing excitation terms as is the case with the other forms of ex- 
citation that are considered in this report. 
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VI. GENERALIZED FEEDLINE COMPUTER CODE 


A versatile computer code has been generated to calculate 
the frequency response of propellant feedlines. Program philosophy 
and logic are discussed in this section, and the mechanics of setting- 
up and execution including specific example problems are presented 
in Appendix H. 

The construction of a computer code to generate the frequency 
response of a propellant feedline is a relatively straightforward task 
when the sequence of feedline components (lines, bellows, PVC joints, 
etc. ) is fixed by a specific design configuration. However, for pre- 
liminary design analyses, it is advantageous to determine the effect 
that a reordering of line components would produce on the overall 
feedline response. For example, various sequences of components 
could materially affect the structural-hydraulic coupling that is 
known to influence the instabilities that were delineated in Section I. 

It can readily be visualized that reprogramming the feedline re- 
sponse equations for numerous changes in the component sequence 
is an extremely expensive, time-consuming task. To alleviate this 
situation, a master computer code was written in FORTRAN IV for 
the CDC 6400 to provide the systems engineer with a design tool for 
obtaining the frequency response of a feedline in which the type, 
number and sequence of basic line components between the propellant 
tank and the turbopump inlet can be specified in a completely arbi- 
trary fashion. The dominant criterion employed in generating this 
code was that the "user" be required to perform a minimal amount 
of manual conditioning of a given problem prior to machine execution. 

Formulation of the code is based on the fact that the four- 
terminal pressure-flow relationship across any line component can 
be represented in matrix form when the perturbation equations are 
transformed into the Laplace domain. That is, with no loss in 
generality. 


■Pi+1~ 


■p r 


= Ri 


-Qi+i. 
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x — 1, 2, . . . . n 


(79) 


where n represents the number of components in the line. In this 
expression, the transformed output pressure and flow for a given 
component are related to the corresponding input quantities through 
a 2 X 2 square matrix, D* , plus a column matrix, Ci , that is 
present only if that particular component is being excited by an 
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external forcing function, K. The parameter, K, may be the 
Laplace transform of an externally imposed acceleration or velo- 
city, or the volume perturbation due to a side branch pulser. In 
general, it is desired to determine the perturbation pressure at 
some point in the line in response to any one of the class of forcing 
functions that K represents. With regard to propellant feedline 
analyses, the important pressure point is located at the entrance to 
the turbopump because dynamic variations in the inducer inlet suc- 
tion head are instrumental in the growth of structural -fluid dynamic 
instabilities. Having established the desired output variable, the 
functional form of the transfer function for the overall line is ob- 
tained by applying Equation (79) to each line component followed by 
sequential matrix substitution to arrive at the following generalized 
transfer equation: 


P = D Q ± K B (80) 

Through experience we have found that matrix D can be stated 
a priori for any feedline composed of n ordered components. 

D = _D n D n _i D n _ 2 .... Dj, (81) 

Matrix .B, however, does not possess such a well-defined property. 
In general, 


B — B L + B g + . . . . + B m (82) 

The value of M m M is strongly dependent on the line configuration as 
is the matrix structure of each Bj . The distinguishing feature of 
each Bj is that it consists of a column matrix of one of the Ci's 
pre-multiplied by one or more _Di’s. Matrix j? in Equation (80) 
contains the desired pressure response to the perturbation, K, and 
matrix Q contains the pressure-flow perturbations at the exit to 
the fuel tank. By assuming that the impedance at the tank exit, Zi , 
is known and that the turbopump-injector -engine combination can be 
lumped into an equivalent impedance, Z t , Equation (80) can be ex- 
panded into its constituent equations to obtain the explicit form of 
the transfer function 

*ji ± jJdgLgl + d a si b n ~ ( d n z i + d i2) ] /oo| 

K (d^Zi + d^) - (d u Z 1 + d 12 )/Z t 

where the bij and d AJ are the elements of matrix B and D, 
respectively. 
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To illustrate the concepts that have been presented, consider 
the feedline shown in Figure 19# It is desired to determine the 
transfer function relating the pressure at the turbopump inlet, P 6 , 
to the line excitation velocity, V,j. Adjacent to each line compo- 
nent is a matrix expression of the type presented in Equation (79). 
Performing the previously indicated matrix substitution yields 


‘ Pe ' 


■ZiQf 


= D 5 Qj, P3 P3 Pi 


Pe/Zt. 


- Ox . 


+ V n (D5 C4 - D5 D4 C3 - D5 D4 P3 


C a ) 


( 84 ) 


The subscripts on pressure and flow rate pertain to specific points in 
the line while the subscripts on matrices _Di and refer to the 

component location in the line. Equation (84) completely defines the 
generalized matrices in Equation (80). An expression similar to 
Equation (83) follows directly from Equation (84). 

The foregoing discussion forms the basis for constructing a 
computer code to determine the frequency response of a feedline con- 
taining an arbitrary number and sequence of components. The com- 
puter code contains a controller program and a separate subroutine 
for each component in which the elements of each D* and Q are 
calculated. The program has the capability of synthesizing a feedline 
with as many as fifteen different types of components. However, the 
program listing that is presented in Appendix H contains eleven 
component subroutines; the four vacancies are available for the 
addition of components during future programs. In addition, the 
source deck contains two subroutines for constructing matrices 13 
and D and one subroutine for calculating the Bessel functions J x 
and J 0 and their ratio for complex arguments. 

Another subroutine calculates the speed of sound in a liquid 
that may or may not contain a homogeneous distribution of bubbles. 

In either case, the radial elasticity effects of the line may be taken 
into account through the Korteweg equation. 

Since substitution of iuo for the Laplace variable n s M implies 
calculations in the complex plane, liberal use has been made of the 
complex variable operations afforded by FORTRAN IV. Complex 
quantities have been subscripted with as many as three indices to 
facilitate bookkeeping within the program. 

In summary, the user merely defines the feedline structure 
as was done in Figure 19 and indexes the components sequentially be- 
ginning with the component that is attached directly to the propellant 
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Station No. Component No. Matrix Pressure Flow Expressions 
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Figure 19. Matrix Representation Of A Feedline In The Laplace Domain 



tank. An expression of the type presented in Equation (80) is then 
assigned to each component. From this point, a minimal amount of 
matrix manipulation is then required to arrive at an expression 
similar to Equation (84), but which reflects the particular configura- 
tion under consideration. The flow of subsequent calculations in the 
program is governed by the following inputs: 

(1) A coded set of integers describing the type of elements 
in the order in which they appear in the line, 

(2) The number of Bi's that are to be summed to arrive at 
matrix and 

(3) The number and type of matrices that are contained in 
each B^i . 

Using the codes in Item (1), additional data, relevant to each compo- 
nent type, are read in and associated with a specific component lo- 
cation in the line. Categorically, this data can include component 
lengths, line radii, friction factors, bubble radii, spring constants, 
damping constants, liquid-vapor mass ratios, etc. The complete 
input package also includes information that pertains to the line in 
general, e. g. , input and terminal impedances, thermodynamic and 
fluid dynamic properties of the propellant, mean flow velocity, the 
elastic and geometric properties of the conduit wall, and the fre- 
quency band over which the response is to be determined. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


VII. 1. Conclusions 


The objective of this study has been to develop an analytical 
model and corresponding computer program to allow for the study 
of disturbances of liquid propellants in typical engine feedline systems. 
This model has been completed and includes (a) the effect of steady 
turbulent mean flow, (b) the influence of distributed compliances, 
such as dissolved ullage gases and flexible walls, (c) the effects of 
local compliances, such as cavitation regions and complex side 
branches, and (d) various factors causing structural-hydraulic 
coupling, such as bends, mounting stiffness, forced changes in line 
length and bellows or PVC joints. The computer program has been 
set up such that the amplitude and phase of the terminal pressure/ 
input excitation is calculated over any desired frequency range for 
an arbitrary assembly of various feedline components. A user f s 
manual has been prepared and is attached as Appendix H of this 
Interim Report. 

Investigation has shown that the predominant effect of turbu- 
lence is to increase the spatial attenuation at low frequencies; at 
high frequencies, the laminar and turbulent frequencies coincide. 

The effect of turbulence also has a tendency to reduce the phase 
velocity at low frequencies; however, this effect can be neglected 
for all feedlines where the parameter r 0 2 U)/v > 10 4 , which includes 
virtually all cases of interest. An additional factor, r c , has been 
added to the laminar propagation operator to account for the turbu- 
lent attenuation contribution. 

For the homogeneous distribution of very small bubbles en- 
trained in the liquid propellant, the net effect over the frequency 
range of interest is to lower the phase velocity of the pure liquid. 

A constant-quality model has been assumed for both single- 
component and two- component, two-phase flows. 

It has been concluded that the effect of the wall compliance 
can be correctly modeled by the classic Korteweg correction to the 
phase velocity; the attenuation for the elastic wall is virtually in- 
distinguishable from the attenuation of a rigid wall. We have also 
concluded that for typical feedline problems, the axial wall stiff- 
ness effect will not permit real wave propagation of the higher order 
modes since the amount of energy which is fed into the higher order 
modes at the line termination will not permit sustained coupling of 
the fluid and the wall for these modes. 
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No other specific conclusions as to the accuracy or appli- 
cability of the feedline model can be made until the completion of 
the experimental study currently in progress at Southwest Research 
Institute. 

VII. 2 Recommendations 


Recommendations for future work include: 

(1) The development of a subroutine to compute the speed 

of sound in two phase, single- component fluids, using the 
equilibrium model described in Section III. 1 of this report. 
This would involve programming a vast amount of thermo- 
dynamic data for each liquid propellant under considera- 
tion. 

(2) Improvement of the model used to compute the bubble 
dynamics of a large local cavitation bubble. The effec- 
tive bubble spring rate (or bubble compliance) is also 
dependent upon the vaporization and condensation rates 
as the bubble undergoes successive perturbations in 
pressure. 

(3) The continuation of the experimental program currently 
in progress. This testing phase of the program will be 
used to evaluate the applicability of the analytical model 
and computer program described in this report. 
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APPENDIX A 

RATIONAL APPROXIMATE MODEL FOR 
DISTRIBUTED PARAMETER SYSTEMS 


A- 1 



The four -terminal representation of a transmission line that 
has been utilized throughout this report is an exact model for the 
zeroeth mode transfer function equations. The value of this repre- 
sentation for frequency response analyses has been proven. In 
many problems, however, the time domain solution is required, and 
the Laplace transform representation constitutes a convenient tech- 
nique for obtaining the space- dependent solution independently of the 
time -dependency. In the final analysis, the Laplace inversion inte- 
gral must be evaluated. When hyperbolic functions are involved in 
the inversion, as they are in the distributed parameter line model, 
an infinite series of time- dependent terms is obtained due to the un- 
bounded number of zeros of the hyperbolic functions. In many cases, 
closure of this series is extremely difficult, and the desired accur- 
acy does not warrant laborious mathematical operations. The ob- 
vious conclusion is that there is a definite need for a simple, accu- 
rate, approximate technique for obtaining the time domain solution 
from the frequency domain representation. 


The rudiments of such a technique were first proposed by 
Oldenberger and Goodson^, and were later refined and put in a 
practical form by Gerlach^> This technique is based upon ex- 
panding the hyperbolic functions into an infinite product of second- 
order polynomials, i. e. , 
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(A. 1) 


sinh r(s) = r(s) 
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(A. 2) 


The coefficients, C cn , » Csn> a.re evaluated by solving for the 

value of Laplace operator, s n at the zeroes of the hyperbolic func- 
tions. These quantities are conveniently catalogued in Figures A-l 
to A-4. To evaluate any one of these coefficients, it is sufficient to 
calculate the damping number of the line, D n . The damping number 
and the desired term number in the expansion completely define the 
four coefficients. It may be shown that for approximation purposes 
the term T(s), in the expansion for the hyperbolic sine, can be 
represented adequately by sL/c 0 . Characteristic impedance can 
be approximated by p 0 c 0 . 


With these expansions, it is now possible to approximate the 
hyperbolic functions to any desired degree of accuracy with expres- 
sions that can easily be (1) inverted into a time solution using 
standard transform pairs or (2) integrated in the time domain to 
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With Axial Damping Number 




A-5 


With Axial Damping Number 
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Figure A-4. Variation of the Approximate Model Parameter 
C sn With Axial Damping Number 




obtain the real time solution. For example, consider a transfer 
function in which a one-term approximation for the hyperbolic 
functions produces 


P(s) 
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Standard inversion pair yields 
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As a second example, consider a transfer function of the type 

P(s) = -Z e (s) tanh T(s) V(s) (A. 5) 

A one -term approximation in the model yields 



Interpreting the Laplace operator as a time derivative, the follow- 
ing equation is obtained: 
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(A. 7) 


With a second equation relating pressure and flow velocity, a stan- 
dard numerical integration easily can be obtained on the computer. 


The accuracy of this model is summarized below: 

(1) One term of the model well approximates the hyper- 
bolic operators up to the first critical frequency. 

(2) Two terms improve the approximation up to the first 
critical point and roughly (not well) approximate the 
hyperbolic operators up beyond the second critical 
frequency. The use of more terms would improve the 
results near the second critical frequency. 

(3) A one-term approximation for the hyperbolic function 
is generally adequate for most engineering problems. 
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SPEED OF SOUND IN A LIQUID CONTAINING 
A HOMOGENEOUS DISTRIBUTION OF BUBBLES 



Single- Component, Two-Phase Mixture in Equilibrium 


The speed of sound in a pure, two-phase substance can be 
expressed by the standard form 


a 8 = go 
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or, in terms of specific volume rather than density. 
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Following the approach of Gouse and Brown®, it can be shown that 
the partial derivative in Equation (B.2) can be expanded such that 
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Since constant pressure and temperature lines are coincident in the 
two-phase region beneath the vapor dome. 
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or in finite difference form 
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The subscripts, f and g, refer to the saturated liquid and vapor 
state, respectively. Substituting Equations (B.4) and (B.5) into 
Equation (B.6) 


(%i 


(B.6) 


Equation (B.6) can be developed into a form that is suitable 
for use with a standard temperature -entropy chart. The quality of 
a two-phase state, x, is defined as 


x = M. v /(M £ +M v ) (B . 7) 

and the specific volume for such a state can be written in terms of 
quality as 
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V = Vf + X (v g - Vf ) 


(B.8) 


Substituting this expression into Equation (B.6), the speed of sound 
can be written as a function of quality 
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In finite difference form, the speed of sound in a single- component, 
two -phase mixture in thermodynamic equilibrium becomes 


a 




(B. 10) 


Single -Component or Two -Component, Two-Phase Mixture With 
Constant Quality 

The above section described the acoustic velocity in an 
equilibrium mixture of a liquid and its vapor phase. In this section, 
the speed of sound is analyzed for a constant quality mixture of a 
liquid and a distinctly different gas, but the results also apply to 
single -component mixtures as long as the vaporization and conden- 
sation rates are low enough that no phase change occurs as a result 
of pressure disturbances (constant quality is assumed at all times). 
In the development of a constant quality model for the acoustic velo- 
city in a liquid-gas mixture, the following assumptions are made: 

(1) The mixture is homogeneous in phase composition. 

(2) The mass of each phase remains constant. 

(3) The gas behaves as a perfect gas with the appropriate 
equation of state and constant specific heats. 

(4) The liquid compressibility is ignored. 

(5) The gas and liquid phases are always at the same 
temperature. 

The subsequent derivation is based on total volume as opposed to the 
specific volume approach that was utilized in the previous section. 
Assuming an adiabatic process. 
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For a two- component mixture. 
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From the definition for density, 

V = M/p (B. 14) 
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The mass ratio, 0, is now defined as the mass of vapor /mass of 
liquid. 


0 = Mg/Mi 

Substituting the mass ratio into Equation (B.15), 


Vi+ V g = 
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Thus, the speed of sound in a constant quality mixture becomes 
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As mentioned in Section III, the isothermal velocity of sound 
might be more appropriate for two- component, two-phase mixtures 
rather than the adiabatic (isentropic) velocity of sound. 

Plesset and Hsieh^ have shown that bubble dynamics are 
governed by an isothermal process at low excitation frequencies and 
an adiabatic process at higher frequencies, providing there is a uni- 
form temperature distribution within the gas bubbles. However, 
this sharp division of thermodynamic behavior is not necessarily 
the case because the heat conduction rate of the gas is considerably 
smaller than that in the liquid. The liquid has a large specific heat 
and thermal conductivity, and behaves as a heat reservoir. Conse- 
quently, in the liquid adjacent to the gas -liquid interface, there are 
no changes in the gas temperature during compression. In the 
center of the bubble away from a substance having a high specific 
heat, the gas follows the adiabatic equation of state. Therefore, the 
overall bubble follows a polytropic process. However, for very 
small bubbles, the heat diffusion into the liquid is so r^pid that 
temperature changes in the bubble can not take place, and the pro- 
cess is essentially isothermal. 


For this reason, the propagation of sound in a liquid with a 
homogeneous distribution of bubbles (basic assumption No. 1) indi- 
cates a speed in agreement with isothermal speed of sound. To 
compute the isothermal speed of sound in the mixture, the acoustic 
velocity of the gas in Equation (B.13) should be changed to the iso- 
thermal value, or 


c g =Vgo Y RT' 

where y = 1.0 rather than 1.4. 


(B.19) 
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APPENDIX C 

COMPLIANCE OF LARGE BUBBLES AND THE 
POLY TROPIC EXPONENT 


C-l 


4 


The "complex" compliance of large cavitation bubbles was 
defined in Section IV. The bubble compliance was shown to be fre- 
quency dependent, especially near the bubble resonant frequency. 
This section describes how the bubble damping and polytropic ex- 
ponent can be estimated for an arbitrary bubble size and excitation 
frequency. 


In the literature, there is considerable disagreement as to 
the theoretical expressions that should be used to evaluate the 
damping constant, 6, which is related to the coefficient of damp- 
ing, b, by 


6 = 


b 

UJm 2 


(C.l) 


The reason for this disagreement is that, while the theory is appli- 
cable to all frequencies, experimental verification has been obtained 
only at the resonant condition. However, it is hypothesized that the 
theory which is summarized here and has been incorporated into the 
computer code is valid over a broad range of frequencies. 

Total damping, due to the presence of a bubble, is attributed 
to the losses originating from three processes: 

(1) Thermal damping resulting from the thermal conduction 
between the gas in the bubble and the surrounding liquid; 

(2) Sound radiation damping caused by energy dispersed by 
radiating spherical sound waves when the bubble is ex- 
cited into volume pulsations; 

(3) Viscous damping from viscous forces at the liquid-gas 
interface. 


The total damping constant is given by 

6 = &th 4 ^rad 4 &vis 
h = bth 4" b r ad 4 bvla 

where 


bth = 


sinh (20i R 0 ) + sin (20!R O ) 1 \ 

cosh(20xR o ) - cos(20iR o ) 0iR o 
20i R 0 sinh (20iR o ) - sin (20xR o ) / 
3(y-l) cosh(20iR o ) - cos(20i,R o ) ) 


X 


t)Pq 

V o 0) 


(C.2) 


(C.3) 
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PlUJ 5 
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bvis 



In these expressions, 


(C.5) 


|i = liquid viscosity 
= liquid density 
C£ = liquid speed of sound 

Y = ratio of specific heats irrespective of the thermo- 
dynamic process indicated by the polytropic expo- 
nent, T), 


Also, the argument 



(C.6) 


is a measure of the rate at which heat is conducted over a distance, 
R 0 , from the bubble center to the liquid-gas interface. The thermal 
diffusivity of the entrained gas is denoted by Dx . 


The correct value of the polytropic exponent, T), depends 
on the bubble size and excitation frequency. The value of the poly- 
tropic exponent has been found to be 


Y [ 1 + 6th 2 1 

Ti 3(Y-1) / smhj[20j L R 2 ) - sin (2gjRo) \ ] 

L + 20 x Ro \ cosh(20 l R o ) - cos( 203 .R o ) ) J 

where 6 th is equal to the term in brackets in Equation (C.3). 
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APPENDIX D 


FORCED CHANGES IN LINE LENGTH 



The POGO phenomenon is a coupled dynamic instability in- 
volving the vehicle structure, propulsion system and propellant 
feed system. As a direct result of this coupling, the structural 
excitation applied to the feedline may result in a relative motion of 
the ends of the line, i. e. , forced changes in line length. The pur- 
pose of this section is to analyze the changes in length from the 
viewpoint of an inhomogeneous , boundary- valued problem and to 
derive the form of the transfer matrices that describe the pressure- 
flow response to this type of excitation. 

The physical problem is shown in Figure D-l. Liquid pro- 
pellant is flowing through a constant area duct of length L, internal 
radius r c , and wall thickness t. The terminal ends of the line are 
acted upon by structural velocities having arbitrary relative phase 
and magnitude. The relative motion due to this excitation is assumed 
to be small in the linearized sense. 
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Figure D-l. Geometry For Forced Changes In Line Length 

The following assumptions and assertions are made with 
regard to the fluid flow: 

(1) The flow is viscous, incompressible and axi symmetric; 

(2) Viscosity remains constant; 

(3) The flow is laminar; this assumption is not a severe 
constraint because the added attenuation due to turbu- 
lence can be adequately taken into account through the 
use of the modified propagation operator described in 
Section II; 



(4) The pressure and the axial and radial components of 
flow consist of a steady- state level plus a perturbation 
value; 

(5) The axial velocity perturbation, u, is much larger 
than the radial perturbation, v; 

(6) Velocity gradients in the axial direction contribute a 
negligible amount to the viscous dissipation forces; 

(7) Spatial derivatives of the density can be neglected since 
the compressibility of the liquid propellant is extremely 
low. 

D' Souza and Oldenberger^ have shown that when these assumptions 
are applied to the Navier - Stoke' s equations in cylindrical coordi- 
nates, the continuity equation and the equation of state, the follow- 
ing linearized differential equations describe the flow field: 


chj dp r d 2 u 

Po sF = " 5 X + Ldpr 


1 dp dv r v r du 

— — *- + — - + — + — = 
X dt dr r dx 



0 


Su I 
dr - 


(D.l) 

(D.2) 


These expressions are easily recognizable as the equation of motion 
in the axial direction and the combined continuity and state equations, 
where x is the fluid bulk modulus of elasticity. In addition, ab- 
sence of an equation of motion in the radial direction implies that 
the pressure is constant across the cross section, i. e. , p = p(x, t). 


The relative motion in the line wall is governed by the wave 
equation 



where y is the axial displacement of a point in the wall and c„ is 
the longitudinal wave speed in the wall. Equation (D.3) is subject to 
the imposed boundary conditions 


dy(o, t) 
dt 


Vi(t) 


and 


dy(L, t) 

dt 


V a(t) 


(D.4) 


The no- slip compatibility condition at the wall-fluid interface is 
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d^t ~' = u ^ x ’ r °’ 


(D.5) 


The solution for the flow field proceeds in the Laplace trans- 
form domain noting that the average pressure perturbation at an 
arbitrary cross section, p(x, t), is equal to the actual perturbation, 
p(x, t). The transform of the equation of motion with zero initial 
conditions yields 


5 s U 1 SU s / 1 dp\ 

S r 2 r dr v \ p 0 s dx / 


(D.6) 


where the upper case letters indicate transformed variables. Since 
P is independent of the radial coordinate, Equation (D.6) can be re- 
duced to an expression involving only one dependent variable by the 
transformation of variables 


V 



dP 

dx 


(D.7) 


Combining Equations (D.7) and (D.6) yields a Bessel type differential 
equation with one finite solution at r = 0. 

V = f(x, s) Jo (5r) (D.8) 

where §= \f- s/ v and f(x, s) is an undefined function which must be 
determined. Therefore, the transformed velocity, U, becomes 

U = f(x, s) Jo (Sr) - ~ (D.9) 

Solution of the wave equation in the wall subject to the prescribed 
boundary conditions yields 


„ - sx ^ n sx 

Y(x, s) = Bi cosh — + B 2 smh — 
C w C w 


(DAO) 


where B x 


Vi(s) - p V 2 (s)-V 1 (s)coshsL/G 1 
and B 2 = T ; “ 

s s smh sL/C w 


Combining the transformed compatibility condition. Equations (D.5) 
and (D.10), produces 


1 dp 
p 0 s dx 


f(x, s) Jo(?r 0 ) 


[ sx 
Bi cosh — + B 2 sinh 
C M 



(D.ll) 
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Substituting this result into Equation (D.9) and averaging over the 
cross section yields 

U = f(x, s) ^^ — - J 0 (§r 0 )] + s[Bx coshg + B s sinh^] (D.12) 


As before, barred quantities indicate average values. If the con- 
tinuity equation is transformed and then averaged across the cross 
section, the result is 


p = __H su 

s dx 


(D.13) 


Note that the dependence on radial velocity vanishes since there is 
no net radial flow. Differentiating this result with respect to x and 
incorporating Equations (D.12) and D. 11) yields a second-order, in- 
homogeneous, ordinary differential equation for the arbitrary func- 
tion f(x, s). After considerable algebraic manipulation, the equa- 
tion takes the following form: 


d 3 f 

dx 2 


- y 3 f 


A.c £\ g , y 

\ C w s yj 0 (§r 0 ) 


g(x, s) 


(D. 14) 


where 


Y 3 


2 


r 2 f, 2 Ji(?r q) 

0 L ' ?r 0 J 0 (Sr 0 ) 


and 


sx sx 

g(x, s) = Bx cosh — + B 2 sinh — 


A particular integral of Equation (D. 14) is easily obtained by the 
method of undetermined coefficients because of the cyclic nature of 
the inhomogeneous function g(x, s). The general solution of Equa- 
tion (D # 14) is 

f(x, s) = (Ai cosh Yx+A 2 sinh Yx) *f [q coshy^ + C 2 sinh~ ] (D.15) 

\ C w C w J 


where Ai and A 2 are undetermined functions of the Laplace variable 
s and Ci and C 2 are related to Bi and B 2 by 


Ci 


(c^~ - y3 ) Jo(§ r o) 


i = 1, 2 
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After simplification, the transformed velocity field. Equation (D.12) 
becomes 


U = (Aj coshYx+ A 2 sinh Yx) J 0 (§r 0 )j + 

s (^-' y3 )/ 

p 2 r f Bi cosh — + B 2 sinh 

s- - Y s 1 ' C " 

'W ) 


(D.16) 


Noting that K = P 0 c 0 , the pressure field is obtained from Equation 

( D. 1 3 ) . 


P_ -P° c ° 


(Aj. YsinhYx + As Ycosh Yx) - J 0 (§r 0 )^ 


2 

S *,2 

s| 7T5- - Y 

C ° / l „ s . , SX s sx 

B, — sinh — + B 2 — cosh — 


(£ -+) 1 ‘ C< 


(D. 17) 


Coefficients A x and A 2 can be eliminated by evaluating the pressure 
and velocity at the extremities of the line, i. e. , 


P(o, s, ) = Pi 
U(o. s. ) = U 4 


P(L, s) = P 3 
U(L, s) = U 2 


(D. 18) 


In the process of making these evaluations, the characteristic impe- 
dance, Z c , can be defined in terms of the propagation operator, Y, 


Z e = 


/■N 2, . 

Po c 0 Y 


Applying boundary conditions (D.18) to Equations (D.17) and (D.16) 
incorporating the expressions for B x and B 2 , and imposing the 
condition that V 2 (s) = G(s) Vj.(s) yields 


Pa' 


coshYL - Z 0 sinhYL 


Pi 

(JL 

ICo 3 
/ a 2 

- Y 2 ) 

\ Vi(s) 

a i 



-i 



( 

- Y 2 ) 


U 2J 


-Z„ sinhYL coshYL 


LUxJ 

Vc„ 8 

__ 


(D.19) 


where 
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Ot p ~ — 


cosher- - coshYL] + 

G w 


\ (G-cosh-^Y v 

\ + sinh ^ - 7T- sinh Y l\ 

I sinh^ V C “ C « Y > 

C w 


OCi = 


P o C c 


~T~ sinh-^Y -YsinhYL] + 

G w Cvf / 


( G - c ° sh ^) 


sinh 


sL 


— ( cosh—r^ - cosh YlJ 


C w ^ 


Cw 


/ 


The assumed dependence of V 2 (s) on Vi(s) is entirely realistic since 
the vehicle structural velocities, which are the forcing functions, are 
generally describable by transfer functions. 

Equation (D. 19) represents the desired result for forced changes 
in line length. The column matrix in Equation (D. 19) represents a 
pressure-flow modification to the square matrix that describes pres- 
sure and flow in the absence of external excitation. It can be shown, 
though not conveniently, that when G(s) = 1 (rigid body excitation) 
Equation (D. 19) degenerates into the form presented in Section V. 
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APPENDIX E 
MOUNTING STIFFNESS 
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Discrete Parameter - Acceleration Technique 


A typical example of a line with a discrete parameter mount- 
ing stiffness is shown in Figure E-la. The elasticity of the sur- 
rounding structure has been simulated by a spring in parallel with 
a viscous damper. A structural acceleration input, ai(t), mea- 
sured with respect to vehicle fixed axes, acts upon the viscoelastic 
support. Due to the filtering characteristics of the support, the 
applied acceleration produces a response, a^(t), of the vertical 
line segment. It was shown in Section V. 1 that the four-terminal 
representation of a vertically accelerated line is 


Ps“ 


coshYL -Z 0 A 1 sinhYL 


'Pa" 


Z 0 sinh YL 


= 




ajfc(s) 

S 


Cb. 


-AZ 0 1 sinh YE coshYL_ 

| 

_Qa_ 

1 

_A(1 -coshYL) 


Obviously, is identical to ai if the support is infinitely rigid and 

aj % approaches zero if the mount is extremely flexible. In the inter- 
mediate range of interest, a^ can be related to ai through the equi- 
valent spring-mass-dashpot analog shown in Figure E-lb. During 
oscillatory motion, the mass of fluid contained in the transverse limbs 
of the feedline contributes an added mass, (mj + ma), effect. The 
inertial loading, due to the fluid mass in the vertical segment of the 
line, is replaced by the equivalent dynamic pressure forces. For a 
support of negligible mass, the differential equation of motion for 
mass M about its equilibrium position is 

Mx = -k(x-y) - b(x-y) + A(p 3 -p 2 ) (E.2) 

Transforming this result into the Laplace domain and noting that 

^{x} = x, ^ty}=Y 

and 

s 2 X = ajj(s), s 2 Y = ai(s) 


the following result is obtained. 

aji(s) = a i(s)[ Ms 2 +sb+k ] + A t p 3(s)-P 2 (s)] [ Ms 2 +sb+k ] 


(E.3) 


The first term in Equation (E.3) reflects the externally imposed 
acceleration, and the second term reflects a passive inertial loading. 
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Figure E-l. Model For A Line With Mounting Stiffness 



A straightforward combination of Equations (E.3) and (E.l) produces, 
after some manipulation. 


■p 3 ‘ 


h u±£ 

1 


’Pa' 


a 



l+e 

1+3 



a i(s) 





j— - gy ! 



" 1 + 3' 


-Q3-J 


1+3' 

1 +e' . 


-Q»- 


- 0 ? - 


where the (i, j = 1, 2) are the elements of the original square 
matrix in Equation (E.l) and 


a 


sb+k 


s(Ms s +sb+k) 

As 

(Ms 2 +sb+k) 


Z e sinh YE 


Z c sinh yL 


a 


■ .(m^w) A(1 - coshYL > 

A 2 s(l-cosh YL) 

M^+sb+k 


All other terms have been defined previously. The input accelera- 
tion, a^s), could be replaced by the equivalent velocity represen- 
tation, s V i (s). However, this replacement would require a minor 
computer code modification since the program currently determines 
the response to acceleration, ai(s). 


Impedance Technique 

Consider a feedline configuration similar to the one in 
Figure E-la with the exception that the viscoelastic support is re- 
placed by an equivalent driving point impedance, Z 6 . Equation 
(E.l) applies to the vertical line segment. In this case, however, 
we choose to replace the coefficient of the third matrix by its equi- 
valent velocity, Vi(s). The support exerts a force, F, on the 
feedline, and the equation of motion for this system that corresponds 
to Equation (E.2) is 

Mxj a = -F + A(p 3 -p s ) (E.5) 

Applying the Laplace transformation to this expression yields 

Ms V^(s) = -F(s) + A[P 3 (s)-P 2 (s)] (E.6) 
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The driving point impedance is, by definition. 


Z,(s) 

Eliminating 

Vi(s) 


F(s) 

Vji(s) 

F(s) between Equations (E.6) and (E.7) produces 

= A(P a -P s ) 

Z s + M s 


(E.7) 


(E.8) 


Equation (E.8) may be introduced into Equation (E.l), and the follow- 
ing result is obtained after rearrangement: 


'P 3 “ 


in 

i 

.Liz.- 

l+a 


"P 3 " 

-Qs- 


L La+ 1+a 

l+a J 


-Qs. 


where 

AZ C2 sinhr 2 

a = — 

Z s + Mg 

A 8 (l-coshr 2 ) 2 
p ' Zg + Mg 


Lu = elements of original square matrix in 
Equation (E.l). 

The obvious difference between Equations (E.4) and (E.9) is 
that the latter does not contain the trailing column vector that re- 
flects the enforced excitation. In the latter formulation the effect 
of the mounting support (structural impedance) is absorbed into the 
square matrix. Both formulations degenerate into the simple sta- 
tionary line representation for an infinitely stiff support. 
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APPENDIX F 
PARALLEL LINES 
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Consider the assemblage of parallel lines shown in Figure 
F-l. Although this configuration is not used frequently in pro- 
pellant feed systems, the pressure-flow relationships for this com- 
ponent have been developed and included in the computer code for 
application in special situations. Assuming that flow losses due to 
elbow resistance are negligible, a standard four-terminal pressure- 
flow equation can be written for the i-th line. For this particular 
problem, however, it is more convenient to express that pressure- 
flow equation as 

Qj.* = Zcj 1 Ai (Pj coth I* - PgCschTi) 

Q2 t = Zcj 1 At (PiCschTi - PgcothTi) (F.l) 

The terminal pressures, Pj, and P 8 , are common to all line seg- 
ments. The continuity equation requires that 

n n 

£ Qli = Qi an d = Q 2 (F.2) 

i=l i=l 

Summing Equations (F.l) and utilizing Equations (F.2) yields 

n n 

Qi = Pi £ Zqj 1 Ai coth Ii - P 2 £ Zoj 1 Ai csch 
i=l i=l 

n n 

Q 2 = Pi £ Z^At csch T t - P 2 Y, Zc"i l Ai coth r t (F.3) 

i=l i=l 

Equations (F .3) can be easily rearranged into the standard form 


p 2 - 


b u 

^12 


■pr 

Qs. 


-t>2l 

t>32. 


-Qu 


(F.4) 


where 
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Figure F-l. Parallel Line Model 



bn = b; 


n 

£ Zcj 1 Ai coth Tj 
i= 1 


£ Z Ct l Aj c s ch Ft 
i= 1 


bi2 


1 

n 

£ Zoj 1 At csch Ft 
i= 1 


bn 


n 


E z «"i l 

i=l 


A t csch Ft - 


Z" 1 At coth Tt 


£ Z 0i l At csch Tt 
i=l 


The line length that appears in the expression for the propagation 
operator is taken to be the developed length measured along the 
pipe centerline. 
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APPENDIX G 


COMPLEX SIDE BRANCH 
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In general, any type of complicated side element may be 
modeled. These side elements can be visualized as local or lumped 
"compliances, " even though they may not be compliance -like ele- 
ments. The model used for an arbitrary complex side element, 
such as a gauge and line connected onto a feedline, may be put in the 
form of a local compliance, as shown in Figure G-l. This element 
is assumed to consist of a side branch line having an inertance, I, 
and resistance, R, and with a capacitance termination, C. The 
line inertance, I, and the line resistance, R, have been pro- 
grammed respectively ass 

I = p 0 L/A b (G.l) 

and 

R ■ 128 ^ (G - 2) 


Here, p 0 = fluid density, L = side branch length, and A b 
is the branch cross sectional area; |a is the fluid viscosity, and d b 
is the side branch diameter. To consider the effects of compressi- 
bility in the side branch, the inertia and resistance effects are neg- 
lected, and 





for gases 


qt 



for liquids 


(G.3) 


where q* is the flow into the branch, V is the volume of the com- 
pressible fluid contained within the side branch, Y is the ratio of 
specific heats for the fluid (if a gas), and n is the liquid bulk 
modulus. 


For the case of a gauge connected to a feedline containing a 
liquid by a branch line partially filled with air, the lumped model is 
obtained by summing the effects of the branch inertance, resistance, 
and capacitance. 

The resultant expressions in the Laplace domain relating the 
pressure and flow upstream and downstream of the side branch 
become 
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(a) PHYSICAL SITUATION (b) MODEL 



CONTINUITY IN LINE 

Q, (s ) - Qj(s) * Q,(s) 

GAUGE LINE ( Assumed simple inertance plus resistance only ) 

Q 3 <s > = Q 4 (s) 

P e ( s > - R* ( s > - IsQ, + RQ, 

CAPACITANCE 

Q 4 (s) * C s P 4 
"COMPLIANCE MODEL" 

Q,(s) - Qj(s) - 1 IC 2+ RC - j] s P, (s) 

L 5 5 J 2751 

Figure G-l.Example Treatment Of More General "Local Compliance" In Feed Line 
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Pa(s) = Px(s) 


Qs(s) - Ql(s) - [ ICs 2 +RCs + 1 ] sP l(s) 

The quantity in the brackets can be visualized as a complex 
"compliance" for this case, and represents a local compliance 
so far as the feedline is concerned* 


(G.4) 


G-4 



APPENDIX H 

USER'S MANUAL FOR GENERALIZED 
FEEDLINE PROGRAM 
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H. 1 Program Organization 

This manual describes the mechanics of the Generalized 
Feedline Computer Code. It is assumed that the user has had pre- 
vious experience in FORTRAN IV coding. 

The FORTRAN IV source deck is composed of a main or 
controller program, individual subroutines for each line component 
and several special purpose subroutines which either perform 
specific mathematical operations or serve as function evaluators. 
The subroutine approach was chosen because of its distinct advan- 
tages in the areas of program debugging, modification and expan- 
sion. Using the Control Data Corporation 6400 system, the source 
deck is compiled each time the program is submitted to the terminal 
for execution, thus eliminating tape handling. The function of each 
routine is presented below in the order in which it will appear in the 
FORTRAN listing which is included in the next section. 

Controller: 


The controller program performs the following functions: 


(1) Read- in of all pertinent data for a particular line 
configuration; 

(2) Execute the calling of subroutines in the proper sequence; 

(3) Evaluate the magnitude of the transfer function at each 
frequency; 

(4) Increment frequency, the independent variable, within 
a prescribed range; 

(5) Perform all output operations. 


Subroutine ONE : Calculates matrix elements for a straight duct 

of finite length and constant cross-sectional 
area having no external excitation. 


Subroutine TWO; Calculates matrix elements for a single cavi- 

tation bubble including the effects of radiation, 
thermal and viscous damping. 

Subroutine THREE: Calculates the matrix elements of a pressure- 

volume compensator (PVC) joint including the 
effects of internal friction. 
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Subroutine FOUR; Establishes the continuity requirements for 

a side branch pulser. 

Subroutine FIVE: Calculates the matrix elements for a straight 

line which undergoes rigid body motion as a 
result of an external velocity excitation. 

Subroutine SIX ; Calculates the matrix elements for several 

lines in parallel all of which have a common 
entrance -exit point. 

Subroutine SEVEN: Calculates the matrix elements for a bellows 

with gas trap liner. Relative motion between 
the ends of the bellows is permitted. 

Subroutine EIGHT : Calculates the matrix elements for a line with 

mounting stiffness having a structural accelera- 
tion applied to the viscoelastic support. 

Subroutine NINE: Calculates the matrix elements for a line 

undergoing forced changes in length. 

Subroutine TEN : Calculates matrix elements for a line with 

mounting stiffness using the impedance method. 

Subroutine ELEVEN ; Calculates the matrix elements for a com- 
plex side branch with laminar flow. 

Subroutine TWELVE 

to FIFTEEN : Available for future expansion. 

Subroutine SPEED : Calculates the speed of sound in a pure liquid 

or a liquid with a homogeneous distribution of 
an ullage gas. A mixture of a liquid and its 
own vapor may also be considered. However, 
in the latter case, the results are not precise 
because the theoretical model does not in- 
clude the thermodynamic effect of phase change 
at the bubble surface. Finally, the speed of 
sound may be corrected for wall elasticity. 

Subroutine DMAT: Calculates the resultant matrix D in 

Equation (81). 

Subroutine BMAT; Calculates the resultant matrix B in 

Equation (82). 
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Subroutine J1ORJ10 : Evaluates Bessel functions Ji(Z) and J 0 (Z) 

for complex Z. Series expansion is used 
for moderate values of the argument; an 
asymptotic approximation is used for large 
arguments. 

H. 2. Program Listing 

This section contains a listing of the program source deck. 
Program control cards, which are described in Section H. 4, are 
applicable to CDC 6000 series computer systems. 

H. 3. Program Input Package 

Due to the generalized nature of this computer program, the 
length of the input package is variable and is directly related to the 
number of components in the feedline model. The first card in the 
package is always an alphanumeric header card on which the user 
can punch pertinent information such as run number, configuration 
number, etc. , to be printed as the first line of program output. 
Fifty-five (55) Holerith fields have been allocated for the informa- 
tion on the header card. The actual data input is a combination of 
integer or floating point information. Integer data conforms to a 
(2413) format, while the floating point data is read in with a (6E12.6) 
format. 


The following data must be supplied in the order and with the 
units shown: 

CARD 1 Header card 


CARD 2 


NELM, dimensionless Designates the number of 

components in the line. 

JBNUM, dimensionless Number of matrices that 

must be summer to arrive 
at matrix B, Equation (82). 


NGAS, dimensionless Indicates presence (NGAS = 1) 

or absence (NGAS = 0) of dis- 
solved gases in propellant. 

NELAST, dimensionless Includes (NELAST = 1) or 

excludes (NELAST = 0) line 
elasticity effects in speed 
of sound calculation. 
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PROGRAM CONTROLCINPUT, OUTPUT, TAPE bO=lNPUT) 

C NELMrNUMBER OF DISCRETE ELEMENTS IN THE FEEDLINE# I .E. # LINES# 

C BELLOWS # PVC JOINTS#ETC. 

C NELM INTEGER CONSTANTS ARE READ IN AND DESCRIBE THE TYPE OF 

C ELEMENTS IN THE ORDER THEY APPEAR IN THE LINE BEGINNING AT THE 

C OUTPUT OF THE FUEL TANK AND TERMINATING WITH THE COMBINED INPUT 

C IMPEDANCE TO THE TIJRBOPUMP# INJECTOR AND ENGINE. FOR EXAMPLE# 

C ITYPE(J)=M, JsELEMENT NUMBER# M=ELEMENT TYPE 

C Mr 1 SIMPLE LINE WITH NO EXTERNAL EXCITATION 

C Mr? CAVITATION BUBBLE 

C Mr 3 PRESSURE VOLUME COMPENSATOR 

C Mr* SIDE BRANCH PULSER 

C Mr 5 RIGID BODY MOTION OF A SIMPLE L I NE # VELOC I T Y EXCITATION 

C Mrb NPAR LINES IN PARALLEL WITH A COMMON INPUT-OUTPUT POINT 

C Mr? BELLOWS WITH RELATIVE MOTION BETWEEN THE ENDS 

C Mr 8 LINE WITH MOUNTING STIFFNESS(OPTION 1) 

C Mrq FORCED CHANGES IN LINE LENGTH 

C M = 10 LINE WITH MOUTING STIFFNESSCOPTION g) 

C M = ll COMPLEX SIDE BRANCH 

C M = 1 2 

C Mr 1 3 

C M = l* 

C Mr 1 5 

C JBNUMsNUMBER OF MATRICES THAT MUST BE SUMMED TO CONSTRUCT MATRIX B 

C B=B(1)+B(2)+ +BCJBNUM) 

C JTERM(J) IS THE NUMBER OF SUBMATRICES IN B(J) 

C K ( J# M) IS AN ARRAY THAT DESCRIBES THE TYPE OF MATRICES THAT ARE TO 

C BE MULTIPLIED TOGETHER TO CONSTRUCT THE COMPONENTS OF MATRIX B. 

C M TAKES ON VALUES FROM 1 THROUGH JTERM(J) 

C FOR EXAMPLE# IF B(2)rDS*D**C3 THEN K(2#l)r5#K(2#2)=*#K(2#3)r3 

C NGAS=0# NO DISSOLVED GASES IN THE PROPELLANT 

C NGASr 1 , DISSOLVED GASES IN PROPELLANT WITH GAS-TO-LIQUID MASS 

C FRACTION#PHI 

C NELASTrO# SPEED OF SOUND IN PROPELLANT IS BASED ON AN INFINITELY 

C RIGID TRANSMISSION LINE 

C NELASTri# SPEED OF SOUND IN PROPELLANT REFLECTS LINE ELASTICITY 

C (WALL MODULUS EWALL AND WALL THICKNESS HWALL) 

C BRCOMPrCOMPLIANCE OF SIDE BRANCH ELEMENT 

C BRL=LENGTH OF SIDE BRANCH ELEMENT IN FEET 

C BDI AM=DI AMETER OF SIDE BRANCH ELEMENT IN INCHES 

C BRR=RESISTANCE OF SIDE BRANCH ELEMENT 

c ★*★*★*****★*★*★★★*★★★**★**★*★**★★ 

COMPLEX ETA#ENUM#DENOM#TRANS#B#D#XI#ARG#RJ#GAMMA#ZC#COSHG#SINHG#DD 
1 # CC # TCOTH, TCSCH # BB# Z1 #Z2, Xl #X2#X3#Xf, QUADRAT, ALPHAP, BETAP, ALPHADP# 
0BETADP#GOFS#SLCW#COSHCW#SINHCW#COEFF, ALPHA 1# ALPHAB 

complex zstruct#alpha#beta 

DIMENSION ITYPE(15)#JTERM(15)#K(15#1B)#EL(15)#RADIUS(15)#RBUB(15)# 
1FPVC (15) # AKBPVC(15)# AKCPVC ( 15) # PARLEN ( 15 # 5 ) # P ARRAD( 15 # 5 ) # FREQ ( 200 ) 
DIMENSION B(B,i),D(B#B)#SIZE(BOO)#SIZFDB(gOO)#ANGLE(800)#FBEL(15)# 
1C0MPLY(15)#AKBEL(1S),BSIGN(10)#DAMPER(15)#SPRINGK(15) 

DIMENSION DD(g#g#15) #CC(B#l#15) #BB(g# 1#15) # R ADSEC ( 20 0 ) # AREA(iS) 
COMMON PI # RADIUS# OMEGA #ENU, ETA# EL# CO #RHO0#THERMK#RBUB#CPCV#P0# VISC 
1#FPVC# AKBPVC# AKCPVC# NPAR, PARRAD, PARLEN# BULKM0D#RH0LIQ#NGAS#GASMW 
COMMON GAMGAS,G#TO#PHI#COMPLY#NELAST,EWALL#HWALL#NELM, JBNUM# JTERM, 
1K#DD#CC#D#B#BB#CP#J# AREA# VMEAN,FBEL# A KBEL#BSIGN# DAMPER# SPRINGK 
COMMON G1#G2#RH0WALL#ZX,ZY 
COMMON BRCOMP#BRL#BRDIAM 
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********************************* 

1 READ inoo 

IF(EOF,hO)S,10 
5 STOP 
10 CONTINUE 

READ 1005, NELM, JBNUM, NGAS, NELAST 
READ 1005, (ITYPE(J), J=1,NELM) 

READ 1005, CJTERMC J), J=l, JBNUM) 

DO 15 J=l, JBNUM 
KOUNT=JTERM(J) 

15 READ 1005, (KCJ,M) ,M=1,K0UNT) 

READ 1020, CBSIGN(J) ,J=1, JBNUM) 

READ 1020, HERTZ Ir DELHZi , FlL I M, DELHZ2 , HERTZF 
READ 1020, SIGN, TERMZ,RHOLIQ, THERMK,PO, TO 
READ 1020, EWALL,HWALL,ZINPUT,CPCV,BULKMOD, PHI 
READ 1020,VISC,GASMW,GAMGAS,CP, VMEAN, BRCOMP 
C LOOP FOR READ-IN OF LINE ELEMENT DATA 

Pl=3. 1*15927 
DO 95 J = 1 , NELM 
INDEX=ITYPE( J) 

GO TO (20, 25, 30, 35, *0, *5, 50, 55, b0,b5, 70, 75, 80, 85, 90), INDEX 
20 READ 1020,EL(J),RADIUS(J) 

AREA(J)=PI*RADIUS(J)**2/lt*« 

GO TO 95 

25 READ 1020,RBUB(J) ,RADIUS(J) 

GO TO 95 

30 READ 1020, FP VC ( J) , AKBP VC ( J) , AKCPVC ( J) 

GO TO 95 
35 GO TO 95 

“0 READ 1020,EL(J),RADIUS(J) 

AREA(J)=PI*RADIUS(J)**2/1**. 

GO TO 95 

*5 READ 1005, NPAR 

READ 1020, (PARLEN(J, I) ,PARRAD( J,I), 1=1, NPAR) 

GO TO 95 

50 READ 1020,FBEL(J),COMPLY(J), AKBEL(J) 

GO TO 95 

55 READ 1020, EL( J ) , R AD I US ( J ) , 0 A MPER ( J) , SPR I NGK ( J ) 

AREA( J)=PI*RADIUS(J)**2/lf *. 

GO TO 95 

bO READ 1020,EL(J),RADIUS(J),Gl,G2,RHOWALL 
AREA(J)=PI*RADIU3(J)**2/lf * . 

GO TO 95 

K5 READ 1020,EL(J),RADIUSCJ),ZX,ZY 
AREA(J)=PI*RADIUS(J) **2/1**. 

GO TO 95 

70 READ 1020, BRL , BRD I AM 
GO TO 95 
75 READ 1000 
GO TO 95 
80 READ 1000 
GO TO 95 
85 READ 1000 
GO TO 95 
90 READ 1000 
95 CONTINUE 

C ********************************* 

6=32.17*0*9 
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ETA=(0.,1.) 

IFREQ=1 

FREQ(IFREQ)=HERTZI 
CALCULATION OF MATRIX ELEMENTS 
SR 0MEGA=2.*PI*FREQ(IFREQ) 

00 1?5 J=1,NELM 
INDEX=ITYPE(J) 

GO TO ( 100 ,105,110,115, 120 , 125, 130 , 135,1*0, 145, 150,155, lbOrlbS, 170 
l)r INDEX 
100 CALL ONE 
GO TO 175 
105 CALL TWO 
GO TO 175 
110 CALL THREE 
GO TO 175 
115 CALL FOUR 
GO TO 175 
120 CALL FIVE 
GO TO 175 
125 CALL SIX 
GO TO 175 
130 CALL SEVEN 
GO TO 175 
135 CALL EIGHT 
GO TO 175 
DO CALL NINE 
GO TO 175 
D5 CALL TEN 
GO TO 175 
150 CALL ELEVEN 
GO TO 175 
155 CALL TWELVE 
GO TO 175 
150 CALL THIRTEN 
GO TO 175 
155 CALL FOURTEN 
GO TO 175 
170 CALL FIVETEN 
175 CONTINUE 

********************************* 
CALL OMAT 
CALL BMAT 

ENUM=B(i,i)*(D(2,i)*ZlNPUT+D(2,2))-B(2,l)*(D(lrl)*ZINPUT+D(lr2)) 

DEN0M=D(2rl)*ZINPUT+D(2r2)-(D(l, 1 ) *Z I NPUT+D ( 1 r 2 ) ) / TERMZ 

trans=sign*enum/denqm 

SIZE( I FREQ) =CABS (TRANS) 

SIZEDB(IFREQ)=20.*AL0G10(SIZE(IFREQ) ) 

ANGLE ( I FREQ ) = AT AN2 ( A I MAG (TRANS) r RE ALC TRANS)) *5 7.23578 
180 ANGLE(IFREQ)=ANGLE(IFREG)-180. 

185 RADSEC(IFREQ)=OMEGA 

IFCFREQ(IFREQ).LT.HERTZF) GO TU ISO 
GO TO 205 
ISO IFREQ=IFREQ+1 

IF(FREQ(IFREQ-1)-F1LIM)13S,13S, 200 
135 FREQ(IFREQ)=FREQ(IFREQ-1)+DELHZ1 

ifinal=ifreg 

GO TO 3S 

200 FREQ(IFREQ)=FREQ(IFREQ-1)+DFLHZ2 
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IFINAL=IFREQ 
GO TO 99 
205 PRINT 1000 
PRINT 1030 
DO 250 1=1, IFINAL 

250 PRINT 10 to,RADSEC(I),FREQ(I),SlZE(I),SIZEDB(I),ANGLE(I) 

GO TO 1 

C ***************★*********:******** 

1000 FORMAT ( 55H1 
1005 FORMAT (2*113) 

1020 FORMAT (bEl2. b) 

10 30 FORMAT (ItHOOMEGA-R AD/SEC, t X, 7HFREQ-HZ, 8 X,5HTRANS, ?X, 8HTRANS-DB, bX 
19HANGLE-DEG) 

10 tO F0RMAT(1X,F9.1,F15.1, FI 3. 3, Fit. 2, F15. 1) 

1050 F0RMAT(lHl,Elt.brbEl5.b) 

END 
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SUBROUTINE ONE 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
A STRAIGHT LINE OF LENGTH L WITHOUT EXTERNAL EXCITATION 
COMPLEX ETA, ENUM , DENOM , TRANS, B, D, XI, ARG,RJ, GAMMA, ZC , COSHG , S I NHG , DD 
l,CC,TCOTH,TCSCH,BB,Zl,Z2, Xl , X2,X3,X*, QUADRAT, ALPHAP , BET AP , ALPHA DP , 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA1, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (IS), JTERMQ5), K(15, 15 ) , EL(15 ) , R ADIUS(IS) , RBUB(IS) , 
1FPVCC15), AKBPVCC15), AKCPVCC15) , PARLENf IS , S) , PARR AD (15 , S) , FREQ ( 200 ) 
DIMENSION B(2,1),DC2,2) ,SIZE(200),SIZEDB(200), ANGLE(200),FBEL(15) , 
J COMPLY (15), AKBELC15), BSIGN(IO) , DAMPER C IS ), SPRI NGK ( 15 ) 

DIMENSION DD(2,2,15),CC(2,l,lS),BB(2,l,i5),RADSEC(200),AREA(l5) 
COMMON PI, RADIUS, OMEGA ,ENU, ETA, F.L, CO, RHOO , THERMK , RBU8 , CPC V,PO, VI SC 
1, FP VC, AXBP VC, AKCPVC,NPAR,PARRAD,PARLE'N,BULKM0D,RH0LIQ,NGAS,GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, V MEAN, FBEL,AKBEL,BS IGN, DAMPER, SPRINGK 
COMMON GI,G2,RH0WALL,ZX,ZY 
CALL SPEED 

REYNOLD=VMEAN*B.*RADIUS( J)/(ENU*12.) 

RT=if .*C2.*ENU*0.0055 *REYNOLDw*0.8S)/(RADIUS(J)**2) 

Z1=C0MEGA*EL( J)*ETA)/CO 
Z2=CSQRT((1.,0.)+RT/C0MEGA*ETA)) 

XI=SQRT(OMEGA/ENU)*CSQRTC-ETA) 

ARG=XI*RADIUS( J)/12. 

CALL J10RJ0(ARG,RJ) 

GAMMA=ETA*0MEGA*EL(J)/(C0*CSQRT(C1. ,0.)-2.*RJ/ARG))+REALCZl*Z2) 

ZC = RH00*C0/(CSQRT((1. , 0 . ) -2 . *R J/ ARG) ) 

C0SHG=(CEXP(GAMMA)+CEXP(-GAMMA))/2. 

SINHG=(CEXP(GAMMA)-CEXPC-GAMMA))/2. 

DD(i,l, J)=COSHG 

DD(1,2, J)=-ZC*SINHG/AREA(J) 

DD(2,1,J)=-AREA(J)*SINHG/ZC 

DDC2,2, J)=COSHG 

RETURN 

END 
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SUBROUTINE TWO 

THIS SUBROUTINE CALCULATES THE COMPLIANCE OF AND TRANSFER MATRIX 
ACROSS A SINGLE BUBBLE 

CPCVsRATIO OF SPECIFIC HEATS FOR THE GAS IN THE BUBBLE 
CPsSPECIFIC HEAT AT CONSTANT PRESSURE FOR THE GAS IN THE BUBBLE 
COMPLEX ETA, ENUM,DENOM, TRANS, B,0, XI, ARG,RJ, GAMMA, ZC,COSHG,SINHG, DO 
1,CC, TC0TH,TCSCH,B8,Z1,Z2, X 1 , X? , X3 , X* , QUADRAT , ALPHAP, BETAP, ALPHADP, 
2BETADP, GOFS, SLCW, CQSHCW , SINHCW, COEFF ,ALPHA1,ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(15),JTERM(lS),KCl5,15),EL(15),RADIUS(15)fRBUB(lS), 
1FPVC(15),AKBPVC(15), AKCPVC ( IS ) , P ARLENC IS , 5 ) , P ARRAD ( 15 , S ) , FREQ ( 200 ) 
DIMENSION B(2,1),D(2,2),SIZE(200),SIZEDB(200),ANGLE(200),FBEL(1S), 
1 COMPLY (IS ),AKBEL( 15) , BS IGN ( 10 ), DAMPER ( 15 ) ,. SPR I NGK ( IS ) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,1S),RADSEC(200),AREA(1S) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CQ,RHQO,THERMK,RBUB,CPCV,PO,VISC 
1,FPVC, AKBP VC, AKCPVC, NPAR,P ARRAD, PARLEN , BULKMOD, RHOLIQ , NGAS, GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NE LAST, EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, V MEAN, F BEL, AK8EL,BS IGN, DAMPER, SPR I NGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
GASC0N=?78.*CP*(1.“1./CPCV) 

GASDENsl'f't.*PO/(GASCON*G*(TO + ^bO.)) 

DTHERMsTHERMK/(GASDEN*CP*G) 

PRO=SORT(OMEGA*3faOO./(2.*DTHERM))*RBUB(J)/12. 

IF (PRO .GT. 35. >20, 10 

10 Tlr(SINH(e.*PR0)+SIN(2.*PR0))/(C0SH(2.*PR0)-C0S(2.*PR0)) 
T2=(SINH(2.*PR0)-SINC2.*PR0))/(C0SH(2 # *PR0)-C0S(2.*RR0)) 
T3=2.*PR0/(3.*(CPCV-1.)) 

Tt=(Tl-l./PR0)/(T3tT2) 

P0LY=CPCV/((1.+T2/T3)*(l.+Tt**2)) 

GO TO 30 
20 POLY=CPCV 

30 BU9V0L = *+.*PI*RBUB(J)**3/(3.*172 8.) 

CALL SPEED 

SPRING = POLY*PO*l‘f‘f ./BUB VOL 
AMASS2=RH00/(f „ *P I *RBUB ( J ) / 1 2 . ) 

IF(PR0.GT.35.)50,*f0 
fO BTHERMsT'+aSPRING/OMEGA 
GO TO AO 

50 BTHERM=3.*(CPCV°1. ) *SPR ING/ ( 2 . *PRO*OMEGA ) 
bO BRAD=AMASS2*RBUB(J)*(OMEGA**2)/(12.*CO) 
8VIS=1?28.*VISC/CPI*RBUB(J)**3) 

BDAMP=BTHERM+BRAD+BVIS 

DD(1,1,J)=(1.,0.) 

DD(1,2, J)s(0.,0.) 

DD(2,l,J)s-ETA*0MEGA/(-AMASS2*0MEGA**2+ETA*BDAMP*0MEGA+SPRING) 

DD(2,2,J)=(l.,n.) 

RETURN 

END 
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SUBROUTINE THREE 

C THIS SUBROUTINE CALCULATES THE TRANSFER MATRIX FOR A PRE8SURE- 

C VOLUME COMPENSATOR (PVC JOINT) 

COMPLEX ETA, ENUMrDENOM, TRANS ,B,D, XI, ARG, RJ, GAMMA, ZCf COSHG, SINHG, DD 
If CC, TCOT H, TCSCH, BB, ZI, Z2, XI, X2,X3,XY, QUADRAT, ALPHAP, BET AP,ALPHADP, 
2BETADP,G OF S,SLCW,COSHCW,SINHCW,COEFF,ALPHAl, ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(lS)f JTERM( 15), K( 15,15 ),EL( 15), RADIUS (15) fRBUB(15), 
IF PVC (15) ,AKBPVC(15) , AKCP VC ( 1 5 ) , P ARLEN ( 15 , 5 ) , P ARR AD ( 1 5 , 5 ) , FREQ ( 20 0 ) 
DIMENSION B(2,1),D(2,2),SIZE(200),SIZEDB(200),ANGLE(200),FBEL(15), 
1C0MPLY(.15),AKBEL(15),BSIGN(10) , DAMPER ( 15 ), SPRINGK ( 15 ) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,15),RADSEC(200),AREA(15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, ELfCO,RHOO,THERMK,RBUB,CPCV,PO, VISC 
1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQf NGAS,GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NEL AST, E WALL, HW ALL, NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VME AN, FBEL,AKBEL, BSIGN, DAMPER, SPRINGK 
COMMON Gl,G2,RH0WALLf ZX,ZY 
FRICT=FPVC(J) 

AKB=AKBPVC(J) 

AKC=AKCPVC(J) 

AKPVC=AKB-AKC 
DD(1,1, J)=FRICT 
DD(1,2, J)=(0.,0.) 

DD(2,1, J)=(0.,0.) 

DD(2,2,J)s(l.,0.) 

CC(1,1, J)=(0.,0.) 

CC(2,1, J)=AKPVC 

RETURN 

END 
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SUBROUTINE FOUR 

C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 

C A SIDE BRANCH PULSER 

COMPLEX ETA,ENUM, DENOM, TRANS, B, D, XI ,ARG,RJ, GAMMA, ZC , COSHG, 8INHG, DD 
1,CC, TCOTH,TCSCH,BB, Zl,Z2, XI, X2,X3,X<+, QUADRAT, ALPH AP , BETAP , ALPHADP , 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA 1,ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERMQ5), K(1S,15), EL( 15 ), RADIUS (IS ) , RBUB(IS), 
1FPVC(15),AKBPVC(15),AKCPVC(15),PARLEN(15,5),PARRAD(15,5),FREQ(200) 
DIMENSION B(2,i),D(2,2),SIZE(200),SIZEDB(300),ANGLE(200),FBEL(15), 
1C0MPLY(15),AKBEL(15),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DD(2,2,15) , CC ( 2 , 1 , 15 ) , BB (2 , 1 , 15 ) , R ADSEC ( 200) , AREA (IS ) 
COMMUN PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOQ,THERMK,RBUB,CPCV,PO, VISC 
1,FPVC, AKBPVC, akcpvc,npar,parrad,parlen-,bulkmod,rholiq,ngas,gasmw 

COMMON GAMGAS,G, TO, PHI, COMPLY, NE LAST, ENALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VMEAN,FBEL, AKBEL,BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
DD(1,1, J)=(1.,0.) 

DD(1,2,J)=(0.,0.) 

DD(2,1,J)=(0.,0.) 

DD(2,2,J)=(1.,0.) 

CC(1,1,J)=(0.,0.) 

CC(3,1, J)=(1.,0.) 

RETURN 

END 
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SUBROUTINE FIVE 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRICES 
FOR THE LONGITUDINAL RIGID BODY MOTION OF A SIMPLE LINE OF LENGTH 
L. THE MOTION IS A RESULT OF VELOCITY EXCITATION. 

COMPLEX ETA,ENUM, DENOM, TRANS, B, D, XI, ARG,RJ, GAMMA , ZC , COSHG, SINHG, DD 
1,CC, TCOTH,TCSCH,BB,Zl, Z2 , X 1 , X2 , X3 , X* , QUADRAT , ALPHAP, BET AP, ALPHADP, 

2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHAl, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPEC15), JTERM ( 15 ), K ( 15,15 ) , EL (15 ) , RAD IUS ( 15 ) , RBUB ( 15), 
IFPVC(15) , AKBPVCC15), AKCPVCC15) ,PARLEN(15, 5) , PARR AD ( 15 , 5) , FREQ (200) 
DIMENSION B(2,I) ,D(2,2) , SIZE (BOO) ,SIZEDB(20n) , ANGLE ( 200) , FBEL ( 15 ) , 
1C0MPLY(15),AKBEL(15),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DD(2,2,15) ,CC(2,1,1S) , BBC 2, 1,15) ,RADSEC(2 00) , AREA (15) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOO,THERMK,RBUB,CPCV,PO, vise 
1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQ,NGAS,GASMW 
COMMON G AMG AS, G, TO, PHI, COMPLY, NEL AST, EW ALL, HW ALL, NELM,JBNUM, JTERM, 
JK,DD,CC,D,B,BB, CP, J, AREA, VMEAN, FBEL, AKBEL,BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CALL SPEED 

REYN0LD=VMEAN*2.*RADIUS(J)/(ENU*12.) 

RT=1YY.*(2.*ENU*0.0055*REYN0L0**0.85)/(RADIUS(J)**2) 

Z1=(0MEGA*EL(J)*ETA)/C0 

Z2=CSQRT((1. ,0.)+RT/(0MEGA*ETA)) 

XI=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS( J)/12. 

CALL J10R JO ( ARG, R J) 

GAMMA=ETA*0MEGA*ELCJ)/(C0*CSQRT((1.,0.)-2.*RJ/ARG))+REAL(Z1*Z2) 

ZC = RH 00 *C 0 /(CSQRT(a. , 0 . )- 2 . *R J/ ARG) ) 

C0SHG=(CEXP( GAMMA )+CEXP (-GAMMA) )/2. 

S I NHG= (C EXP ( GAMMA )-CEXP (-GAMMA ))/2. 

DDC1,1, J)=C0SHG 

DD(1,2, J)=-ZC*SINHG/AREA(J) 

DD(2,1,J)=-AREA(J)*SINHG/ZC 
DD(2,2, J)=COSHG 
CCC1,1, J)=ZC*SINHG 
CC(2,1,J)=AREA(J)*((1.,U.) -COSHG) 

RETURN 

END 
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SUBROUTINE SIX 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
AN ARBITRARY NUMBER OF LINES IN PARALLEL ALL OF WHICH EMANATE FROM 
AND CONVERGE ON COMMON INPUT-OUTPUT POINTS. 

COMPLEX ETA,ENUM,DENOM,TRANSfBfD,Xl,ARG,RJ,GAMMA,ZCfCOSHG,SINHG,DD 
If CC, TCOTH,TCSCH,BB, Zl,Z2,Xl, X2,X3, X4, QUADRAT, ALPHAP, BETAP, ALPHADP, 
2BETADP, GOFSr SLCWf COSHCW, SINHCWr COEFFr ALPHAlf ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(15),JTERM(15),K(15,lS),EL(i5)fRADIUS(15),RBUB(15), 
1FPVCC1S), AKBPVC(lS),AKCPVCCiS),PARLENUS,S),PARRAD(lS,5),FREQ(200) 
DIMENSION BC2, 1), D (2 , 2 ) , SI ZE (200 ) , SIZEDB (500 ) , ANGLE ( 200 ) , FBEL ( IS) , 
1C0MPLY(1S),AKBEL(1S),BSIGN(10),DAMPER(15),SPRINGK(1S) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,15),RADSEC(200),AREA(15) 

COMMON p i, r adius, omega, enu, eta, el, co,rhoo,thermk,rbub,cpcv,pq, vise 

If FPVCf AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLlQ,NGAS,GASMW 
COMMON GAMGAS , G, TO fPHI, COMPLY, NELA3T,EW ALL, HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP,J, AREA, VME AN, FBEL, A KBEL,BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
DO S IsifNPAR 
CALL SPEED 

C0=C0/SQRT(1.+BULKM0D*2.*PARRAD( J, I)/(14t.*EWALL*HWALL)) 

SQFT = PI*PARRAD( J, I)**2/m. 

XI=SQRT(OMEGA/ENU)*CSORTC-ETA) 

APG=XI*PARRAD(J, I)/I2. 

CALL J10RJ0(ARG,RJ) 

GAMMA = ETA*OMEGA*PAR|_EN(Jf I)/(C0*CSQRT(C1. f o. )-2.*RJ/ARG)) 
ZC=RHOO*CO/(CSQRT((1. , 0. )-2 . *R J/ ARG) ) 
C0SHG=(CEXP(GAMMA)+CEXP(-GAMMA))/2. 

S I NHG= (CEXP( GAMMA )-CEXP( -GAMMA ))/2. 
TC0TH=TC0TH+SQFT*C0SHG/(SINHG*ZC) 

5 TCSCH=TCSCH-fSQFT/(ZC*SINHG) 

00(1,1, J)=TCOTH/TCSCH 

DDC1, 2, J)=-l e /TCSCH 

DD(2,1,J)=TCSCH-TC0TH**2/TCSCH 

00(2,2, J)sTCOTH/TCSCH 

RETURN 

END 
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SUBROUTINE SEVEN 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX FOR 
A BELLOWS WITH A LUMPED COMPLIANCE REPRESENTING TRAPPED GAS UNDER 
THE LINER PLUS A PERIODIC VOLUMETRIC CHANGE DUE TO RELATIVE MOTION 
OF THE ENDS OF THE BELLOWS 

COMPLEX ETA, ENUM,DENOM, TRANS, B,D, XI, ARG,RJ, GAMMA, ZC,COSHG,SINHG,DD 
1, CC,TCOTH,TCSCH,BB,Zl,Z2, XI , X2,X3,Xf , QUADRAT, ALPHAP,BETAP, ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA1, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(lS),JTERMCiS),K(15,15),ELCl5),RADlUSClS),RBUbCl5), 
1FPVC(15),AKBPVC(15),AKCPVC(15),PARLEN(15,S),PARRAD(15,5),FREQ(200) 
DIMENSION B(2,l),D(2,2),SIZE(200),SIZEDB(2n0),ANGLE(200),FBELClS), 
IC0MPLY(1S),AKBELC15),BSIGNC10),DAMPERC1S),SPRINGK(15) 

DIMENSION DD(2,2,l5),CCC2,l,i5),BB(2,l,15),RADSEC(200), ARE A (15) 
COMMON P I, RADIUS, OMEGA, ENU r ETA, EL, CO , RHOO,THERMK, RBUB, CPC V,PO, VISC 
1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLlQf NGAS,GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VME AN, FBEL,AKBEL,BSI GN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
FRICT=FBELCJ) 

C=COMPLY(J) 

AKB=AKBEL( J) 

DDC1, I, J)=FRICT 
DDC1,2, J)=C0.,0.) 

DDC2,1, J)=-C*ETA*OMEGA 
DD(2,2, J)=(1.,0.) 

CC(1,1, J)=(0.,0.) 

CC (2 , 1 , J)sAKB 

RETURN 

END 
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SUBROUTINE EIGHT 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRICES 
FOR THE AXIAL MOTION OF A VERTICAL LINE WITH MOUNTING STIFFNESS 
ON THE TWO ADJACENT HORIZONTAL LINES. THE TWO STIFFNESS TERMS ARE 
LUMPED INTO ONE EFFECTIVE SPRING IN PARALLEL WITH ONE EFFECTIVE 
VISCOUS DAMPER. THE FORCING FUNCTION IS AN ACCELERATION APPLIED 
TO THE SPRING-DAMPER SUPPORT PLUS A PASSIVE INERTIAL LOADING 
COMPLEX ETA, ENUM,DENOM, TRANS, B,D, XI, ARG,RJ, GAMMA ,ZC, COSHG, SI NHG,DD 
1,CC, TC0TH,TCSCH,BB,Z1,Z2,X1»X2, X3, Xt, QUADRAT, ALPHAP, BETAP, ALPHADP, 
«?BETADP,GOFS,SLCW,COSHCW,SINHCW,COEFF, ALPHA 1, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERM ( IS ) , K C IS , 15) ,EL(15) , R ADIUS( IS ) , RBUB ( 15 ) , 
1FPVCC1S) , AKBPVCQ5), AKCPVCC15) , P ARLEN ( IS , 5 ) , P ARR AD( 1 5 , 5 ) , FREQ (200) 
DIMENSION B(2,l) ,D(2,2),SIZE(2n0),SIZEDB(200) , ANGLE (200) ,FBEL(15), 
1C0MPLY(15),AKBEL(15),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DO (2, 2, 15), CC (2,1, 15), BB (2, 1,15) , RADSEC (200 ) , AREA (15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, C0,RHO0,THERMK, RBUB, CPCV,PO, VISC 
1 , FP VC , AKBP VC , AKCP VC , NPAR ,PARRAD, PARLEN, BULKMOD, RHOL IQ , NGAS , G ASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NELAST,EWALL,HWALL,NELM,J8NUM, JTERM, 
1K,DD,CC,D,B,BB,CP, J, ARE A , VME AN , F BEL , AKBEL , BS I GN , DAMPER , SPR INGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CALL SPEED 

REYN0LD=VMEAN*2.*RADIUS(J)/(ENU*12.) 

RT = l‘f4-.*(2.*ENU*0.0055*REYN0LD**0.85)/(RADIUS(J)**2) 

Z1=(OMEGA*EL(J)*ETA)/CO 

72=CSQRT((1. ,0.)+RT/(0MEGA*ETA)) 

Xl=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/12. 

CALL J10RJ0(ARG,RJ) 

C,AMMA = ETA*OMEGA*EL( J)/(CO*CSQRT( (1. ,0. )-2.*RJ/ARG))+REAL(ZJ.*Z2) 
ZC=RH00*C0/(CSQRT((1.,0.)-2.*RJ/ARG)) 
C0SHG=(CEXP(GAMMA)+CEXP(»GAMMA))/2. 
SlNHG=(CEXP(GAMMA)-CEXP(-GAMMA))/2. 

AMASS=RHOO*(AREA( J-1)*EL( J-1)*AREA( J+1)*EL( J+l) ) 
QUADRAT=-AMASS*OMEGA**e+ETA*OMEGA*DAMPER( J)+SPRINGK( J) 
ALPHAP=(ETA*OMEGA*DAMPER( J)+SPRINGK( J))*ZC*S1NHG/ (QUADRA T*ETA*OMEG 
1A) 

BETAPsAREA( J) *ETA*OMEGA*ZC*SlNHG/QUADRAT 

ALPHADP=(ETA*OMEGA*DAMPER(J)tSPRINGK(J) ) * ARE A ( J) * ( ( 1 . , 0 . ) -COSHG) / ( 
1ETA*0MEGA*QUADRAT) 

RETADP=(AREA(J)**2)*ETA*0MEGA*((1.,0.)-C0SHG)/QUADRAT 
DD(1,1, J)=(C0SHG+BETAP)/(1.+BETAP) 

00(1 ,2, J)=-ZC*SINHG/(AREA(J)*(1.+BETAP)) 

00(2,1, J)=-AREA(J)*SINHG/ZC+BETADP*((1. ,0.) -COSHG) /(1.+ BETAP) 

00(2, 2, J)=C0SHG-f-8ET ADP*ZC*SINHG/ ( AREA ( J ) * ( 1 . 4-BET AP) ) 

CC(1 ,1, J)=ALPHAP/(1.+BETAP) 

CC(2, 1, J)=ALPHADP/(1.+BETAP) 

RETURN 

END 
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SUBROUTINE NINE 

C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX 

C FOR A LINE WITH FORCED CHANGES IN LINE LENGTH 

COMPLEX ETA,ENUM,DENOM,TRANS, B,D, XI, ARG, RJ, GAMMA , ZC , COSHG, SINHG, DD 
If CC,TCOTH, TCSCH, 88, Zl, Z2, Xl,X2,X3, XH, QUADRAT, ALPHAP, BE TAP, ALPHADP, 
2BETADP, GOFSf SLCWfCOSHCW, SINHCWfCOEFF, ALPHAlf ALPHAS 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(IS), JTERM ( 15 ) , K ( IS , IS ) , EL( IS ) , RADI US (15 ) , RBUB (15 ) , 
1FPVC(15), AKBPVC(15), AKCP VC ( 1 5 ) , P ARLEN ( 15 , 5 ) , P ARR AD ( 1 5 , 5 ) , FREQ ( 200 ) 
DIMENSION B(2,1),D(2,2),SIZEC200),SIZEDBC200), ANGLE! 200 ), FBEL ( IS) , 
lCOMPLY(lS) ,AKBELaS),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,I,15),BBC2,1,15),RADSEC(200),AREA(15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO,THERMK, RBUB, CPCV,PO, VISC 
1,FPVC, AKBPVC,AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLlQ,NGAS,GASMW 
COMMON GAMGAS,G, TO, PH I, COMPLY, NEL AST, E WALL, H WALL, NELM,JBNUM, JTERM, 
IK, DD,CC,D,B,BB, CP, J, AREA, VMEAN,FBEL,AKBEL,BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
G0FS=Gl+G2*ETA 

CW=SQRT((EWALL*G)/(RH0WALL*12.)) 

CALL SPEED 

REYN0LD=VMEAN*2,*RADIUS(J)/(ENU*12.) 

RT = l‘f‘f.*(2.*ENU*0.D055*REYN0LD**0.85)/(RADIUS(J)**2) 

Z1=(OMEGA*EL(J)*ETA)/CO 

Z2=CSQRT((1. ,0.)+RT/C0MEGA*ETA)) 

XI=SGRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/12. 

CALL J10RJ0(ARG,RJ) 

GAMMA=ETA*OMEGA*EL(J)/(C0*CSQRT( Cl. ,0.)-2.*RJ/ARG))+REALCZl*Z2) 
ZC=RH00*C0/(CSQRT((1.,0.)-2.*RJ/ARG)) 

COSHG=(CEXP(GAMMA)+CEXP (-GAMMA) )/2. 
SINHG=(CEXPCGAMMA)-CEXP(-GAMMA))/2. 

DD(1,1, J)=COSHG 

DD(1,2, J)=-ZC*SINHG/AREA(J) 

DD(2,1, J)=-AREA(J)*SINHG/ZC 
DD(2,2,J)=C0SHG 
SLCW=ETA*OMEGA*EL(J)/CW 
COSHCW=(CEXP(SLCW)+CEXP(-SLCW) )/2. 
SINHCW=(CEXP(SLCW)-CEXP(-SLCW))/2. 

C0EFF=(-0MEGA**2/C0**2-( GAMMA/EL (J))**2)/(-0MEGA**2/CW**2-( GAMMA/E 
1L(J))**2) 

ALPHAJ. = (RH00*CQ**2/(ETA*0MEGA))*( (ETA*OMEGA*SINHCW/CW-GAMMA*SINHG/ 
1EL( J))+(GOFS-COSHCW)*(COSHCW-COSHG)*ETA*OMEGA/(CW*SINHCW)) 
ALPHA2=-((C0SHCW-C0SHG)+(G0FS-C0SHCW)*(SINHCW-ETA*0MEGA*EL(J)*SINH 
1G/(CW*GAMMA))/8INHCW) 

CC(1,1, J)=C0EFF*ALPHA1 

CC (2 , J , J)=C0EFF*ALPHA2*AREA( J) 

RETURN 

END 
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SUBROUTINE TEN 

C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER MATRIX 

C FOR THE AXIAL MOTION OF' A VERTICAL LINE WITH MOUNTING STIFFNESS. 

C THE STIFFNESS IS REPRESENTED BY THE IMPEDANCE PRESENTED TO THE 

C LINE BY THE VEHICLE STRUCTURE. 

C THIS SUBROUTINE CAN BE USED ONLY WHEN ANOTHER SOURCE OF EXCITATION 

C IS PRESENT IN THE FEEDLINE, I.E. PULSER OR VERTICAL LINE VELOCITY 

C EXCITATION. THE FUNCTIONAL FURM OF THE STRUCTURAL IMPEDANCE IS 

C THAT OF AN EQUIVALENT SPRING AND DASHPOT IN PARALLEL. 

C ZX AND ZY ARE THE REAL AND IMAGINARY PARTS OF ZSTKUCT. 

COMPLEX ETA,ENUM,DENOM,TRANS,B,D,XI,ARG,RJ,GAMMA,ZC,COSHG,SINHG,DD 
1,CC, TC0TH,TC8CH,BB, Zl,Z2, Xl ,X2,X3,Xf , QUADRAT, ALPHAP , BETAP, ALPHADP , 
2BETADP,GOFS,SLCW,COSHCW,8INHCW,COEFF, ALPHA1 , ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPEU5),JTERMC15),K(1S,15),EL(15),RADIUSU5),RBUB(15), 
1FPVCC1S), AKBPVCC15),AKCPVC(15),PARLENUS,5),PARRAD(1S,S), FREQ (200) 
DIMENSION B(2,1),D(2,2),SIZE(20O),SIZEDB(2OO),ANGLE(2OO),FBEL(15), 
1C0MPLYC15), AKBEL ( 15 ),BSIGN(10), DAMPER (IS) ,SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,15),RADSEC(200),AREA(15) 
COMMON PI, RADIUS, OMEGA, ENU, ETA, EL, CO, RH00,THERMK,RBU8,CPCV,P0, VISC 

1 ,FPVC, AKBPVC, akcpvc,npar,parrad,parlen,bulkmod,rholiq,ngas,gasmw 

COMMON GAMGAS,G, TO, PHI, COMPLY, NELAST,EwALL,HWALL,NELM,JBNUM,JTERM, 
IK, DD,CC,D,B,BB, CP, J, AREA, V ME AN, FBEL, AKBEL, BSIGN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
CALL SPEED 

PEYN0LD=VMEAN*2.*RADIUS(J)/(ENU*12.) 

RT = lt'+.*(2.*ENU*0.0055*REYNOLD**0.85)/(RADIUS(J)**2) 

71=(0MEGA*EL(J)*ETA)/C0 

Z2=CSQRT((1.,0.)+RT/(0MEGA*ETA)) 

XI=SQRT(OMEGA/ENU)*CSQRT(-ETA) 

ARG=XI*RADIUS(J)/12. 

CALL J 10R JO ( ARG , R J) 

GAMMA=ETA*OMEGA*EL( J)/(CO*CSQRT( (1. ,0. )-2.*RJ/ARGI ) +REAL C Z 1*Z2 ) 
ZC=RHOO*CO/(CSQRTC(I. ,0. )-2.*RJ/ARG) ) 
C0SHG=(CEXP(GAMMA)+CEXPC-GAMMA))/2. 

SINHG=CCEXP(GAMMA)-CEXP (-GAMMA ))/2. 

AMASS=RHOO*C AREAC J- 1 ) *EL C J-l ) + ARE A C J+l ) *EL C J+l ) ) 

zstruct=zx»et a*zy/omega 

ALPHAS AREAC J)*ZC*SINHG/(ZSTRUCT*AMASS*ETA*OMEGA) 

BETAs(i.-C0SHG)*C AREAC J ) **2 ) / C ZSTRUC T+ AM ASS*E T A*OMEG A ) 

DOCl, 1, J)sCCOSHG+ALPHA)/Cl.+ALPHA) 

DO Cl, 2, J)=-ZC*SINHG/CAREAC J) * Cl . + ALPHA ) ) 

DDC2, 1, J)=-AREAC J) *SI NHG/ZCtBETA* C 1 . -COSHG) / C 1 . t ALPHA ) 

DDC2,2, J)=COSHG+ZC*SINHG*BETA/CAREAC J)*Cl.+ALPHA)) 

RETURN 

END 
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SUBROUTINE; ELEVEN 

THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE TRANSFER 
MATRIX FOR A COMPLEX SIDE BRANCH. THIS MODEL ASSUMES A 
LAMINAR FLOW IN THE SIDE BRANCH. 

COMPLEX ETA,ENUM, DENOM, TRANS, B,D, XI, ARG,RJ, GAMMA, ZC,COSHG, SI NHG,DD 
1,CC, TCOTH,TCSCH, BB,Z1,Z2, Xl , X2, X3 , X* , QUADRAT, ALPHAP, BETAP, ALPHADP, 
2BETADP,G0FS,SLCW,C0SHCW,SINHCW,C0EFF, ALPHA1, ALPHA2 
COMPLEX ZSTRUCTf ALPHA, BETA 

DIMENSION I TYPE (15), JTERM ( 15 ) , K ( 15 , IS ) , EL ( 15 ) , R AD I US( 15 ) ,RBUB(1S) , 
1FPVCC15) , AKBPVC(15),AKCPVC(15),PARLEN(15,5) , PARR AD (15, 5) , FREQ (200) 
DIMENSION B(2,1),D(2,2),SIZE(200),SIZEDB(200),ANGLE(200),FBEL(1S), 
1C0MPLY(15), AKBEL(1S),BSIGN(10),DAMPER(15) ,SPRINGK(15) 

DIMENSION DD(2,2,15) ,CC(2, 1, 15) , BBC 2, 1,15) ,RADSEC(200) , AREA(15) 
COMMON PI,RADIUS,OMEGA,ENU,ETA,EL,CO,RHOO,THERMK,RBUB,CPCV,PO, VISC 
1,FPVC, AKBPVC, AKCPVC,NRAR,PARRAD,PARLEN f BULKMOO, RHOLIQ, NGAS, GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NE LAST, EWALL,HWALL,NELM,JBNUM, JTERM, 
1K,DD,CC,D,B,BB, CP, J, AREA, VMEAN,FBEL,AKBEL,BSIGN, DAMPER, SPRINGK 

COMMON Gl,GE#RHOWALLr ZX,ZY 
COMMON BRC0MP,BRL,8RDIAM 
COMPLEX S 
S=ETA*0MEGA 
10 CALL SPEED 

AREAB=(PI/t.)*(BRDIAM/12.)**2 

BRI=RHOO*BRL/AREAB 

BRRsl28.*VISC*BRL/(PI*(BRDIAM/12.)**f.) 

DD(1,1, J)=(1.,0.) 

DD(1,2,J)=(0.,0.) 

DD(2,1,J)=-S*BRC0MP/(BRI*BRC0MP*3**2+BRR*BRC0MP*S +1.) 
DD(2,2,J)s(l.,0.) 

RETURN 

END 
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SUBROUTINE SPEED 
THIS subroutine calculates THE 

FEEDLINE. THIS FLUID MAY BE A 
A LIQUID AND AN ULLAGE GAS 
COMPLEX ETA,ENUM,DENOM, TRANS, B 
1,CC, TC0TH,TC8CH, BB, Z1,Z2,X1,X2 
2BtTADP,GOFSf SLCW, COSHCWr SINHCW 
COMPLEX ZSTRUCT, ALPHA, BETA 
DIMENSION ITYPE(15), JTERM(15), 

1FPVCUS) , AKBPVCU5), AKCPVC(IS) 

DIMENSION B(2,l) ,0(2,2) , SIZEC2 
1 COMPLY C15) , AKBELC15) ,BSIGN(10) 

DIMENSION DD(2,2,15),CC(2,1,15 
COMMON PI, RADIUS, OMEGA, ENU, ETA 
1, FPVC, AKBPVC,AKCPVC,NPAR,PARRA 
COMMON GAMGAS,G, TO, PHI, COMPLY, 

1K,DD,CC,D,8,BB,CP, j,area,vmean 
COMMON G1,G2,RH0WALL,ZX,ZY 
CLIQUID=SQRT(BULKMOD/RHOLIQ) 

IFCNGAS.EQ.l)i0,20 
10 RHOGAS = l‘f*f.*PO*GASMW/(15H5.*(TO + ‘fbO.)*G) 

CGAS = SQRT(GAMGAS*G*15fS.*(TO + ‘fbO.)/GASMW) 

E1=RH0GAS+PHI*RH0LIQ 

E2=RhOGAS/(RHOLIO*CLIQUID**2*E1)+(PHI*RHOLIQ)/(E1*RHOGAS*CGAS**2) 

E3=(1.+PHI)*RH0GAS*RH0LIQ 

C0=SQRT(E1/(E2*E3)) 

RHOO=E3/EI 
ENU = V I SC/RHOO 
GO TO 30 
20 CO=CLIQUID 
RHOO=RHOLIQ 
ENU=V I SC/RHOO 
30 IF(NELAST.EQ.l)f0,50 

*0 C0=C0/SQRT(l.+BULKM0D*2.*(RADIUS(J)+HWALL)/(lYf . *EWALL*HWALL ) ) 

50 RETURN 
END 


SPEED OF SOUND FOR THE FLUID IN THE 
LIQUID OR A HOMOGENEOUS MIXTURE OF 

, D, XI, ARG,RJ, GAMMA, ZC, COSHG, SINHG,DD 
,X3,Xt, QUADRAT, ALPHAP,BETAP,AL PH ADP, 
,COEFF, ALPHA1,ALPHA2 

K(I5,15),EL(I5),RADIUS(I5),RBUB(15), 
,PARLEN(15,5),PARRAD(15,S),FREQ(200) 
00),SIZEDBC2no),ANGLE(200),FBEL(15), 
, DAMPER ( 15 ),SPRINGK( 15) 

),BB(2, 1, 15) ,RADSEC(200), AREA (15) 
,EL,CO,RHOO, THERMK,RBUB,CPCV,PO, VISC 

d,parlen,bulkmod,rholiq,ngas,gasmw 

NELAST,EWALL,HWALL,NELM, JBNUM, JTERM, 
,FBEL,AKBEL,BSIGN,DAMPER,SPRINGK 
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SUBROUTINE DMAT 

C THIS SUBROUTINE CALCULATES MATRIX D IN THE EXPRESSION P=DQ+VB 

COMPLEX ETA,ENUM, DENOM, TRANS, B,D, XI , ARG,RJ, GAMMA, ZC , COSHG , SINHG, DD 
CCfTCOTHrTCSCH f BBfZlf 22 f XI# XPrXBrXf, QUADRAT, ALPHAP,BETAPrALPHADPf 
2BETADP,GOFS,SLCW,COSHCW,SINHCW,COEFF,ALPHA1,ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION I TYPE (15), JTERM(IS) , K ( 15 , i 5 ) , EL ( 15 ) , RADIUS ( 15 ) , RBUB ( 15 ) , 
1FPVC(15),AKBPVC(15), AKCPVC ( 15) , PARLEN(15 , 5) , P ARRAD (1 5 , 5 ) , FREQ ( 200 ) 
DIMENSION B(2,1),DC2,2) , S I ZE ( 200 ) , SI ZEDB ( 200 ) , ANGLE ( 200 ) ,FBEL( 15) , 
1COMPLY(15),AKBELC15),BSIGN(10),DAMPER(15),SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,1,15),BB(2,1,15),RADSEC(200), AREA (15) 
COMMON P I, RADIUS, OMEGA, ENU, ETA, EL, CO, RHOO,THERMK,RBUB, CPC V,PO, VISC 
1,FPVC, AKBP VC, AKCPVC, NPAR, PARR AD, PARLEN,BULKMOD, RHOL IQ »NGA3,GASMW 
COMMON GAMGAS,G, TO, PHI, COMPLY, NE LAST, EWALL,HWALL,NELM,JBNUM,JTERM, 
1K,DD,CC,D,B,BB,CP,J, AREA, VME AN, FBEL,AKBEL,BS I GN, DAMPER, SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
DO 10 M=1 , 2 
DO 10 N=1 , 2 
10 D(M,N)=DD(M,N,1) 

IF(NELM-1)S0,20,30 
20 RETURN 
30 DO <f0 J=2,NELM 

X1=DD(1,1, J)*D(1,1)+DD(1,2, J)*D(2,1) 

X2=DD(1,1, J)*D (1,2) +DD( 1,2, J) *0(2,2) 

X3=DD(2,1, J)*D(1,1)+DD(2,2, J)*D(2,1) 

XY=DD(2,1, J)*D(1,2)+DD(2,2, J)*D(2,2) 

0(1,1)=X1 
D(1,2)=X2 
D(2,1)=X3 
fO D(2,2) = X*f 
RETURN 
50 PRINT bO 
CALL EXIT 

bO FORMAT (28H0ERR0R IN NELM SPECIFICATION) 

END 
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SUBROUTINE BMAT 

C THIS SUBROUTINE CALCULATES MATRIX B IN THE EXPRESSION P=DQ+VB 

COMPLEX ETA, ENUM,DENOM, TRANS, B,D, XI , ARG,RJ, GAMMA, ZC, COSHG, SINHG, DD 
l,CC,TCOTH,TCSCH,BBf Zl,Z2f XI, X2^X3,XH, QUADRAT, ALPHAP,BETAPfALPHADPf 
2BETADPfG0FS,SLCW,C0SHCW,SINHCW,C0EFFf ALPHA 1, ALPHA2 
COMPLEX ZSTRUCT, ALPHA, BETA 

DIMENSION ITYPE(IS), JTERM( 15) ,K( 15,15), EL( 15), RADIUS (15), RBUB(15), 
1FPVC(15), AKBPVC(15),AKCPVC(15),PARLEN(15,5),PARRAD(1S,5),FREQ(200) 
DIMENSION B(2,1),D(2,2),SIZE(200),SIZEDB(200), ANGLE (200), FBELC15), 
1C0MPLY(1S),AKBEL(15), BSIGN(IO) ,DAMPER(15) , SPRINGK(15) 

DIMENSION DD(2,2,15),CC(2,l,15),BB(2,l,15),RADSEC(200),AREA(15) 
COMMON PI, RADIUS, OMEGA, ENU,ETA,EL, CO, RHOO,THERMK,RBUB,CPCV,PO, VI SC 
1,FPVC, AKBPVC, AKCPVC,NPAR,PARRAD,PARLEN,BULKMOD,RHOLIQr NGAS,GASMW 
COMMON GAMGAS,G, TO , PHI , COMPLY , NEL AST, E WALL , HW ALL, NELM, JBNUM, JTERM, 
1K,DD,CC,D,B,BB,CP, J, AREA, VME AN, FBEL, AKBF.L, BSIGN, DAMPER ,SPRINGK 
COMMON G1,G2,RH0WALL,ZX,ZY 
DO 100 1=1, JBNUM 
KOUNT=JTERM( I ) 

L=K ( I , KOUNT) 

DO 100 J=l, KOUNT 
I F ( J-l ) 5 , 5 , 15 
5 DO 10 M=1 , 2 
10 BB(M,i, I)=CC(M,1,L) 

GO TO 100 

15 L=K(I,K0UNT-J+1) 

Xl=DD(l,l,L)*BB(l,l,I)+DD(l,2,L)*BB(2,l,I) 

X2=DD(2,1,L)*BB(1,1,I)+DD(2,2,L)*BB(2,1,I) 

BB(1,JL, I)aXl 
BB(2 ,1,1 )=X2 
100 CONTINUE 

B(1,1)=BSIGN(1)*BB( 1,1,1) 

B(2,1)=BSIGN(1)*BB(2,1,1) 

110 DU 200 M=1 , 2 

DO 200 1=2, JBNUM 

200 0(M,1)=B(M,1)+BSIGN(I)*BB(M,1, I) 

210 RETURN 
END 
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SUBROUTINE J10RJ0(Z,RJ) 

C CALCULATION OF JO(Z) AND J1(Z) Z-COMPLEX 

COMPLEX Z, J5 , JO,TERMO,TERMl,Zl,Z2,PO,QO,Pi,Ql,PHO,PHl,FZl,EZ2,RJ 
X=REAL(Z) 

Y=AIMAG(Z) 

RsCABSCZJ 

IF(R-18. ) 100, 100, 110 

100 TERMl=Z/2. 

Jl=Z/2. 

J0=(1.,0.) 

TERM0 = (1. ,D.) 

A = 1 . 

AM=15.+R 

101 TERM0=TERM0*(“(Z/2.)**2)/A**2 
JO=JO+TERMO 

TERMl=TERMl*(-CZ/2.)**2)/(CA+l.)*A) 

Jl=Jl+TERMl 

A=A+1. 

IF(A-AM)101,101,115 

110 IF(X)111, 112,112 

111 Zl=-Z 

GO TO 113 

112 Z1=Z 

113 PI=3.1fl5S2b 
Z2=8.*Z1 

I F(CABS(Z2)-5000. 5120,120,121 
12 0 P0=(1. ,0. )-9.5/Z2**2+3b?5./(8.*Z2**'0 
Q0=-l./Z2+37.5/Z2**3-59S3S./(8.*72**5) 

Pl=(l. ,0.)+?.5/Z2**2-f 725./(8.*Z2**'0 
Ql=3./Z2-52.S/Z2**3+72?b5./(8.*Z2**5) 

GO TO 122 

121 P0=(1. ,0. )-t.5/Z2**2 
Q0=-l./Z2 

Pl=(l.,0.)+?.5/Z2**2 

Ql=3./Z2 

122 PH0=Z1-PI/^. 

PH1=Z1-.?5*PI 

FZl=2./(PI*7l) 

FZ2=CSQRT(FZ15 
AZ1 = A I MAG (Zl) 

IFCABS(AZl)-50.)llb,llb,117 
lib J0=FZ2*(P0*CC0S(PH0)-Q0*CSINCPH0)) 

J1=FZ2*CP1*CC0S(PM1)-Q1*CSIN(PH1)) 

IFCXDllf ,115,115 
111 J1=-J1 
115 RJ=J1/J0 
GO TO 118 
11? SIN=Y/ABS(Y) 

RJ=(SIN*(0. ,1* )*P1+Q1)/(P0-SIN*(0. ,1.)*Q05 
119 RETURN 
END 



CARD 3 


ITYPE(J), Jr 1, A sequence of dimensionless 

NELM dimensionless numbers that describes the 

type of elements in the line 
in the order in which they 
appear beginning with the 
element attached to the fuel 
tank outlet. See the first 
page of program listing for 
an expanded definition. 


CARD 4 JTERM (J), Jr 1, The number of submatrices 

JBNUM dimensionless in B (J). For example, in 

Equation (84), JBNUM =3, 
JTERM (1) r 2, JTERM (2) = 3 
and JTERM (3) =4. 


CARD 5 K(J, M), dimensionless 

6 


4 + BJNUM 


An array that fixes the order 
of matrix multiplication in 
B(J). For example, in Equa- 
tion (84), B(3) =^D 4 D 3 C 2 . 
Then K(3,l) = 5, K(3,2) = 4, 
K(3,3) = 3 and K(3,4) = 2. 


CARD B SIGN (J), dimensionless Floating point array that de- 


(5 + JBNUM) 

CARD HERTZI, Hz 

(6 + JBNUM) 

DELHZ1, Hz 

F1LIM, Hz 

DELHZ2, Hz 
HERTZ F, Hz 


fines the algebraic sign pre- 
ceding each B(J). Again in 
Equation (84), BSIGN(l) = 

+ 1.0, BSIGN(2) — 1.0 and 
BSIGN(3) = -1.0. 

The initial frequency where 
response calculations begin. 

Frequency step size for fre- 
quencies between HERTZI 
and F1LIM. 

Frequency break-point which 
defines the boundary between 
regions of different grid size. 

Frequency step size for fre- 
quencies beyond F1LIM. 

The final frequency to be 
evaluated. 
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CARD 

(7 + JBNUM) 


SIGN, dimensionless 


Algebraic sign preceding 
the parameter K in Equa- 
tion (80). 


TERMZ, lb- sec/ft 5 Terminal impedance of the 

line which reflects the com- 
bined effect of the turbopump- 
injector- engine combination. 
In its present form, TERMZ 
is a scalar resistance. 

RHOLIQ, lb-sec 2 /ft 4 Liquid propellant density. 

THERMK, Btu/(hr-ft-°R) Thermal conductivity of the 

gas inside large bubbles. 

Mean propellant pressure. 

Mean propellant temperature. 

Wall modulus of elasticity. 

Wall thickness. 

Scalar flow impedance at the 
propellant tank outlet. 

CPCV, dimensionless The ratio of specific heats 

for the gas in the cavitation 
bubble at the propellant 
temperature. 

BULKMOD, lb /ft 2 Bulk modulus of the pro- 

pellant. 

PHI, dimensionless Vapor-liquid mass ratio for 

the bulk propellant flow. 

CARD VISC, lb- sec/ft 2 Propellant dynamic viscosity. 

(9 + JBNUM) 

GASMW, dimensionless Molecular weight of dissolved 

vapor or ullage gas for use in 
speed of sound calculation. 

GAMGAS, dimensionless Ratio of specific heats for the 

gas whose molecular weight 
is GASMW. 


P0, psia 
TO, °R 

CARD EWALL, lb /in. 2 

(8 + JBNUM) 

HWALL, in. 

Z INPUT, lb -sec /ft 5 
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CP, BtuAb-R Specific heat at constant 

pressure for the gas in the 
cavitation bubble. 

VMEAN, ft/ sec Mean propellant feed velocity. 

BRCOMP, ft 6 / lb Compliance of the complex 

side branch. 

The sequence and number of the remaining input cards cannot 
be stated a priori since they depend on the component structure of the 
feedline. However, the input for each line component appears se- 
quentially in the same order as the components in the line beginning 
with the component attached to the fuel outlet. Listed below are the 


inputs for each class of component. 
Simple line: EL, ft 

RADIUS, in. 


Cavitation bubble: RBUB, in. 

RADIUS, in. 

Pressure- volume FPVC, 

compensator: dimensionless 

AKBPVC, ft 2 
AKCPVC, ft 2 

Side branch pulser: 

Line with velocity EL, ft 

excitation: RADIUS, in. 

Parallel lines: NPAR, 

dimensionless 
PARLEN, ft 
PARRAD, in. 

Bellows with relative FBEL, 
motion: dimensionless 

COMPLY, ft 6 / lb 
AKBEL, ft 3 


Line length. 

Line radius as measured to 
the inner surface. 

Bubble radius. 

Internal radius of the line 
surrounding the bubble. 

Friction factor expressed as 
the fraction of the inlet pres- 
sure remaining at the outlet. 

Volume change constants 
for compensator. 

No input required. 

Same as simple line. 

Same as simple line. 

The number of parallel 
lines. 

Line length. 

Line radius. 

Bellows friction factor, same 
definition as FPVC 
Compliance of bellows. 

Volume change constant. If there 
is no relative motion, AKBEL =0. 
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Line with mounting 

EL, ft 

Line length. 

stiffness (discrete 

parameter): 

RADIUS, in. 

Line radius. 


DAMPER, 

Coefficient of viscous 


lb- sec/ft 

damping. 


SP RIN GK , lb / f t 

Spring constant. 

Forced changes in 

EL, ft 

Line length. 

line length: 


RADIUS, in. 

Line radius. 


G1 , G2 

Real and imaginary parts 

-■ 


of G(s). 


RHOWALL, 
lb /in 3 

Wall density. 

Line with mounting 

EL, ft 

Line length. 

stiffness (impedance 



technique): 

RADIUS, in. 

Line radius. 


ZX, ZY 

Components of structural 
impedance. 

Complex side branch: 

BRL, ft 

Branch line length. 


BRDIAM, in. 

Branch radius. 


H. 4. Deck Structure for a Normal Compilation and Multiple 

Execution on the CPC 6400 Digital Computer 

Deck structure for multiple production runs is shown on the 
following page. 

The anticipated run time, core memory and output linage 
information indicated on the first two control cards are for illustra- 
tion purposes only. Typically, the combined central processor/ 
print processor time for compilation, execution and listing of one 
run is approximately 25 seconds. In the case of production runs, 
time-line charges can be improved by inserting a NOLIST card in 
front of the PROGRAM CONTROL card which will delete the 
source program listing from the output. 
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H. 5. Computer Symbol List 

In order to facilitate program changes by the user, a corre- 
spondence list of the symbols used in the program is given below. 
The quantities on the left are the FORTRAN alphanumeric symbols 
in the program; on the right is the corresponding item from the 
analysis. Dummy variables, which are frequently used to facili- 
tate programming of long equations, have been omitted. 


TERMZ 

-- 

z t 

EL 


L 

RHOLIQ -- p £ 

(Mass density of propellant; 

RBUB 


Ro 

liquid phase) 


PARLEN 


L p 

THERMK 

-- 

k 

PARRAD 


r P 

PO 

-- 

Po 

FBEL 

-- 

F( P V Z 2 ,G) BEL 

TO 


T 0 

COMPLY 


C 

EWALL 

-- 

E t 

AKBEL 

-- 

K B, BEL 

HWALL 

-- 

h 

FPVC 

-- 

F(PV Z 2 , G) pvc 

Z INPUT 

-- 

Zi 

akbpvc 

-- 

k b, pvc 

CPCV 

-- 

c„/c v 

ACKPVC 

-- 

K C, PVC 

BULKMOD 


K 

DAMPER 

-- 

b 

PHI 

-- 

0 

SPRINGK 

-- 

k 

vise 

-- 

n 

OMEGA 

-- 

UU 

GASMW 


MW 

ENU 

-- 

V 

GAMGAS 

-- 

Y 

REYNOLD 


n r 

CP 


Cp 

RT 


Rt 

VMEAN 

-- 

Vo 

ETA 


(+i) 

RADIUS 



r 0 

XI 

_ _ 

? 
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JIORJO 

-- 

Jl/Jo 

GAMMA 


r 

ZC 

-- 

Z« 

COSHG 

-- 

cosh r 

SINHG 

-- 

sinh r 

AREA 


A 

AMASS 

-- 

m* + m3 

ALPHAP 

-- 

a' 

BETAP 

-- 

p' 

ALP HA DP 

-- 

a" 

BETADP 

-- 

P" 

CW 


c w 

GOFS 


G(s) 

RHOWALL 


Pt 

SLCW 


sL/ C w 

COSHCW 


cosh ( sL/ C„) 

SINHCW 


sinh ( sL/ C w ) 

ALPHA1 


°a 

ALPHA2 

-- 

a 2 

C LIQUID 


c 6 


RHOGAS 


Pg 

CGAS 

-- 

c 6 

RHOO 

-- 

Po 

ENU 

-- 

V 

PRO 

-- 

0iR o 

DTHERM 

-- 

Dl 

B THERM 


b th 

BVIS 


b VISC 

BRAD 


b RAD 

POLY 

-- 

ri 

VBUB 

-- 

V BUB 

SPRING 

-- 

■n p o/v BUB 

ZSTRUCT 

-- 

Z s = Z x - iZy/ UJ 

AREAB 


A b 

BRI 

-- 

I 

BRCOMP 


C 

BRR 


R 

BRL 


u 

BRDIAM 

_ _ 

d b 
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H. 6. Example Problems 


This section contains five example problems which will ac- 
quaint the user with the mechanics of program setup and execution. 
The development of each problem follows the same general format: 
statement of the physical problem, generation of the system matrix 
equation, presentation of critical input data and a computer listing 
and graphical representation of the configuration response. In all 
cases, only a portion of the output is shown since the remaining 
calculations follow the same format through termination of execution. 

Problem No. 1: 


The first feedline example, shown in Figure H-l, consists of 
a rigid length of line, a sidebranch pulser and a cavitation bubble in 
series with the fuel tank and a scalar terminal impedance, R. The 
perturbation pressure at Station 1 is assumed to be negligible. This 
assumption is for illustration purposes only, since the computer pro- 
gram accepts an arbitrary scalar input impedance at Station 1. The 
physical distance between Stations 2 and 4 is assumed to be small 
compared to the total line length, L. It is required that the perturba- 
tion pressure at Station 4 be determined for oscillatory pulser in- 
puts, Qd . The first step is to set up the station numbers, component 
numbers and component matrix equations in the Laplace domain as 
indicated in the illustration. Note that the subscripts on pressure 
and flow perturbations correspond to a particular station number, 
whereas the subscripts on matrices Eh and _Ci refer to a specific 
component number. Beginning with the matrix equation for compo- 
nent number 3, successive substitution of the remaining expressions 
yields the overall line equation: 


- p 4 - 


- 0 - 


= p3 p2 a 


_P 4 /R. 


_Qi_ 


+ Q d (s) (+D 3 C 3 ) 


(H.l) 


The fundamental operations involved in the derivation of this system 
equation are common to the setup of all problems, and the user is 
required to perform these steps in order to execute the program. 

At this point, the user would make the following general ob- 
servations regarding the program input: 

(1) NELM = 3 The number of components in the line. 
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Station No. Component No. Component Matrix Equation 
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Figure H-l. Line Model For Example Problem No. 



(2) B = D 3 


Therefore, JBNUM = 1 

JTERM(l) = 2 (two sub- 

matrices in Bj. = B ) 

K(l,l) = 3 
K(l,2) = 2 


(3) ITYPE (1) = 1 
ITYPE (2) = 4 
ITYPE (3) = 2 


ITYPE(J) indicates the type of compo- 
nent corresponding to component 
number J. Refer to comment cards on 
first page program listing. 


(4) SIGN = +1.0 The algebraic sign preceding the exci- 

tation, Qa , in Equation (H. 1). 


(5) BSIGN (J) =1.0, J = JBUM = 1. BSIGN(J) is the algebraic 
sign preceding the Jth term in B = Bj + Bg+ .... The 
usefulness of this input will be demonstrated later in a 
more complex example. 


In this example, we will assume that the propellant does not contain 
any dissolved gases, and the calculated speed of sound in the pro- 
pellant will reflect wall elasticity effects. 

Using typical values for the physical parameters associated 
with the line and the propellant, a complete input package would take 
the form shown in the following listing. At this point, the reader 
should thoroughly review the setup instructions listed at the begin- 
ning of the PROGRAM CONTROL. These instructions pertain to 
program inputs between statement numbers 10 and 15, and they indi- 
cate the flow of calculations within the program. Note that, when a 
pulser is involved in the feedline, no input is required for its des- 
cription (see statement number 35). 


The computer inputs and outputs for this problem are shown 
on the following page. The column labeled TRANS is the magnitude 
of Pi/Qa. Due to the flexibility of the computer code, the transfer 
function is not normalized. The reader may wish to normalize the 
output, and, in this event, a suitable normalization factor would be 
the characteristic line impedance, Z c = P 0 c 0 . The response magni- 
tude for this configuration is presented in Figure H-2. 
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INPUT 


iFEFOLINE CHECKOUT RUN N0.+ 
3 10 1 


+. bOO00Q+02 +.500000+00 

t. 220000+01 


1 + 2 
2 

3 2 

+.100000+01 
+.100000+01 
+.100000+01 
+.300nQO+08 
+ . + 0 8 0 1 ! 0 - 0 5 
+. 3000U0+02 
+ . 1 2 0 n U 0 + 01 


+.500000+00 
+ . >f 51000 + 05 
+.bb0000-0l 
+. ndooou+un 
+. + 00000+01 
+ . + 0 0 o U 0 + 0 1 


+. 000000+00 

+.000000+00 


+.+b000D-02 

+.1+0000+01 

+. 22 + 000+00 


+.b00000+02 

+.3+7000+02 

+.188000+08 

+.500000+02 


-.288000+03 

+.000000+00 


OUTPUT 


FEEDLINE CHECK0UT RUN Nfc'.4 


0MEGA-RAD/SEC FREQ-HZ 

6.3 1.0 

9.4 1.5 

12.6 2.0 

15.7 2.5 

18.8 3.0 

22.0 3.5 

25.1 4.0 

28.3 4,5 

31.4 5.0 

34.6 5.5 

37.7 6.0 

40.8 6,5 

44.0 7.0 

47.1 7.5 

50.3 8.0 

53.4 e.5 

56.5 9.0 

59.7 9.5 

62.8 10.0 

66.0 10.5 

69.1 11.0 

72.3 11.5 

75.4 12.0 

78.5 12.5 

81.7 13.0 

84.8 13.5 

88.0 14.0 

91.1 14.5 

94.2 15.0 

97.4 15,5 

100.5 16.0 

103.7 16,5 

106.8 17.0 

110,0 17.5 


TRANS 

trans-dr 

1 198.086 

141.77 

1810.423 

150.03 

2440.009 

156.00 

3093,494 

160.74 

3778,293 

164.74 

4502.879 

168.25 

5277, 145 

171.42 

6112,872 

174.36 

7 024.328 

177.14 

8029.074 

179.82 

9149,044 

182.43 

10411 .997 

185.01 

11853.479 

187.61 

13519.412 

1 90.24 

15469.358 

192.93 

17780.115 

195.72 

20548,103 

198.61 

23885.499 

201.62 

27696.005 

204,72 

32596.330 

207,84 

3772 6.364 

210.76 

42452.590 

213.12 

45336.366 

214.44 

45244.860 

214.40 

42503.523 

213.15 

38465.238 

211.15 

34282.647 

208.85 

30502.971 

206.51 

27270.078 

204,27 

24556.224 

202.17 

22263.1 24 

200.23 

20369.256 

198.44 

18744.376 

196.77 

17351.841 

195.23 
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Problem No. 2; 


In this example, the feedline configuration consists of a large 
cavitation bubble situated in the line midway between the fuel tank 
outlet and the terminal impedance, R. As shown in Figure H-3, the 
entire line is acted upon by a structural velocity, V^(s), For analysis 
purposes, it is assumed that bubble acceleration effects are negligible, 
and that the bubble diameter is small compared to the total line length, 
La + L 3 . Proceeding as in Problem No. 1 with the assumption of zero 
tank exit pressure perturbation, the following system equation is 
obtained: 


■ P4 ' 


‘ 0 ' 


= D3 Da Dj, 


.P4/R- 


-Qi. 


-Vi( s) 


Ds D3 Cq + C3 


(H.2) 


Note that in this and the preceding example, the number of line com- 
ponents is the same. However, the present configuration results in 
a more complex functional representation of matrix 13: 

B = D3 Da Cg + O3 — Bi + Ba (H. 3 ) 


The following deductions can be made. 


(l) 

NELM = 3 

(2) 

JBNUM = 2 

(3) 

JTERM(l) = 3 


JTERM (2) = 1 

(4) 

K (1, 1) = 3 


K (1,2) = 2 


K (1,3) = 1 


K (2, 1) = 3 

(5) 

ITYPE (1) = 5 


ITYPE (2) = 2 


ITYPE (3) = 5 


The number of elements in the line. 

The number of Bi's that must be 
summed to form B. 

The number of matrices that comprise 
Bj. andBa, respectively. 


Two-dimensional integer array desig- 
nating the ordered subscripts of the 
matrices that comprise Bj. and Bg . 


Integers describing the type of compo- 
nent located in the Jth component 
position. 
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Station No. Component No. Component Matrix Equations 
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Figure H-3. Line Model For Example Problem No„ 



(6) B SIGN ( 1 ) = + 1 . 0 
BSIGN (2) = +1.0 


The algebraic sign preceding the 
submatrices, Bi. 


(7) SIGN = -1.0 The algebraic sign preceding the 

excitation, V^. 

The complete input package, followed by the program output is 
shown on the following page. The transfer function, which is the 
ratio of P 4 and V^, may be normalized with respect to the pro- 
duct of the acoustic impedance, Z 0 , and the line cross-sectional 
area. A, as shown in the response plot of Figure H-4. Note that, _ 
in this example, both line segments have the same cross-sectional 
area; however, this is not a prerequisite. 
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INPUT 


1FEEPLTNE PROGRAM CHECKOUT RUN NO. 5 


3 

5 

J 

3 

3 

+ . 
+ . 
tm 

m 

+ . 
+. 
+. 
+ . 
+ . 


0 

5 


c 

t 

1 

5 1 


mooon+ui 
inonoo+oi 
mono o+oi 

300000+08 
+ 08000-05 
150000+02 
120 00 0 + U 1 
150000+02 

OUTPUT 


,100000+01 
, sooooo + no 
, + 81000 + 05 
b b 0 00 0-01 
000000+00 
+ 00000+01 
, + 00000 + 01 


+. +00000+01 


, bonooo + oe 
, 220000+01 
000000+00 
000000+00 


+.500000+00 +. bOOOOO+02 
+.+b0000-02 +.3+7000+02 
+.1+0000+01 +.153000+08 
+.22+000+00 +.500000+02 


FEEDLINE PR0GRAM CHECK0UT RUN N0,5 


A-RAD/SEC 

FREQ-HZ 

TRANS 


TRANS-DB 

6,3 

1.0 

416,804 


120.65 

9,4 

1.5 

627,164 


128.82 

12,6 

2.0 

840,211 


134.67 

15,7 

2.5 

1056,945 


139.26 

18.8 

3,0 

1278,440 


143.07 

22,0 

3.5 

1505,862 


146.34 

25,1 

4.0 

1740,510 


149.24 

28,3 

4.5 

1983,841 


151.86 

31,4 

5.0 

2237.523 


154.26 

34,6 

5,5 

2503.479 


156.51 

37,7 

6,0 

2783.964 


158.63 

40,8 

6,5 

3081,641 


160.66 

44,0 

7.0 

3399,698 


162.63 

47,1 

7,5 

3741,988 


164.55 

50,3 

8.0 

4113,216 


166.44 

53,4 

8,5 

4519.181 


168.32 

56.5 

9.0 

4967,098 


170.21 

59,7 

9.5 

5466,009 


172.13 

62.8 

10.0 

6027,303 


174.08 

66,0 

10,5 

6665,330 


176.09 

69.1 

11,0 

7398,010 


178.18 

72,3 

11. 5 

8247,116 


180.35 

75,4 

12.0 

9237,282 


182.62 

78,5 

12,5 

10391,344 


184.97 

81,7 

13.0 

11716.354 


187.37 

84,8 

13,5 

13169,506 


189.71 

88,0 

14,0 

14593.687 


191.77 

91,1 

14,5 

15652.858 


193.17 

94,2 

15,0 

15909,707 


193.49 

97,4 

15,5 

15158.146 


192.53 

100,5 

16,0 

13655,134 


190.44 

103,7 

16.5 

11869.681 


187.63 

106,8 

17.0 

10141.404 


184.49 

110,0 

17,5 

8612.917 


181.22 


.238000+03 

. 000000+00 
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Figure H-4. Frequency Response Of The Feedline In Example 

Problem No. 2 
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Problem No. 3: 


Setup procedures for the feedline example shown in Figure 
H-5 parallel quite closely the procedures used in Problem No. 1. 
The basic difference is in the matrix representation of the center 
component, which reflects a structural acceleration input to the 
line through a viscoelastic mounting element. The functional form 
of the system equation, 


-ai(s) D 3 C 2 (H.4) 

is virtually identical to Problem No. 1. Despite this apparent re- 
dundancy, this problem affords the opportunity to incorporate a pre- 
viously unused subroutine. Subroutine EIGHT, into the analysis. 

Since most vehicles employ ^lightly damped structures, a 
rather large amplification factor, Q, was assumed in order to size 
the spring constant and damping coefficient. By definition. 



The mass term is taken to be the combined fluid mass contained in 
the horizontal limbs of the feedline. For the line dimensions shown 
in Figure H-5 and Q = 100, k and b were taken to be 88 X 10 4 
lb/ft and 45.2 lb-sec/ft, respectively. 

For this problem, 

(1) ITYPE (1) = 1; 

ITYPE (2) = 8; 

ITYPE (3) = 1; 

(2) SIGN = -1.0 

Input values for items (1), (2) and (5) in Problem No. 1 apply also to 
this example. The response magnitude, normalized with respect to 
BOX density and a characteristic length, is shown in Figure H-6. 


P 4 

p 4 /rJ 


= D3 Dg Dj 


LOi J 
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INPUT 


1LINE WITH MOUNTING STIFFNESS SWRI PROJECT 05-2889 

3 1. 0 1 

1 8 1 
2 

3 2 


+.luonoo+ui 
+. 100000+01 
- . 1 0 0 0 0 0 + 0 1 
+ .3000 0 0 + IJ 8 
+. + 080 00- 05 
+. 150000+02 
+ . 1 U 0 n 0 n + 0 2 
+. 150000+02 


+.500000+00 
+.+blOOO+OS 
+. bbOOOO-Ol 
+.000000+00 
+. + 00000+01 
+. + 00000+01 
+. + 00000+01 


+. b00000+02 
+.220000+01 
+.000000+00 
+.000000+00 

+. +52000+02 


+.500000+00 

+. 000000+00 

+.000000+00 

+.000000+00 

+ . 880000 + 01) 


+. b00000+02 
+.3+7000+02 
+.199000+08 
+.500000+02 


OUTPUT 


LINE WITH MOUNTING STIFFNESS SWRI PROJECT 02-2889 


OMEGA-RAD/SEC 

FREO-HZ 

TRANS 

TRANS-OB 

h.3 

1.0 

22.1211 

h 1 . 9 9 

9.+ 

1.5 

22,22b 

b 2 . 0 3 

1-2 . b 

P.n 

22.393 

b 2 . 1 7 

15.7 

2.5 

22 . bl 9 

b 2 . 3 8 

18.8 

3.0 

22.90? 

b 2 . b 3 

22.0 

3.5 

2 3 , 2 b 0 

b2 . 9 + 

25.1 

+ .0 

2 3 . b 8 3 

b 3 . 3 0 

2R.3 

+ .5 

2+.181 

b 3 . ? 1 

31. + 

5.0 

2 + . 7 b 0 

b + . 1 8 

3 + # b 

5.5 

25. + 28 

b+.?2 

37.7 

b.O 

2b. 19b 

b 5 . 3 1 

+ 0.8 

b.5 

27.07+ 

b 5 . 9 7 

+ + .0 

7.0 

28.0?b 

b b . 7 0 

+ 7.1 

7.5 

29.218 

b 7 . 5 0 

50.3 

8.0 

30.519 

b 8 . 3 7 

59.+ 

8.5 

32.000 

b 9 . 3 1 

5b. 5 

9.0 

3 3 . b 8 b 

70.3 + 

59.7 

9.5 

35 . b05 

71. + 5 

bp. 8 

in.o 

37,78b 

72. b + 

bb.O 

10.5 

+ 0 . 2b 0 

73.91 

b 9 , 1 

11.0 

+3.050 

75.25 

72.3 

11.5 

+ b . 1 b 5 

7b. b + 

?5. + 

12.0 

+9.985 

78,07 

7R.5 

12.5 

53.230 

79, +9 

81.7 

1 3 . 0 

5b. 939 

80.8 + 

00 

* 

oc 

13.5 

bO.+Ob 

82.02 

88.0 

1 + .0 

H3.2++ 

82.9 + 

91.1 

1 + • 5 

b + . 9 9 b 

83. +9 

9 + . 2 

15.0 

b5. 328 

83.59 

97.+ 

15.5 

b+ . 1 72 

83.23 

100.5 

lb.O 

bl .759 

82. + 7 

103.7 

lb. 5 

58.5+5 

81. + 0 

10b. 8 

17.0 

5+ . 9+ 3 

80.13 

lin.o 

17.5 

51,308 

78. ?b 


.298000+03 

. 000000+00 
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Figure H-6. Frequency Response Of The Feedline In Example 

Problem No. 3 
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Problem No. 4: 


In this fourth example, shown in Figure H-7, an externally 
excited segment of line is coupled at its extremities to the remainder 
of the line through two flexible bellows. This problem was dis- 
cussed briefly in Section VI. 

Initially, the user would set up the physical problem exactly 
as shown in Figure H-7. Note that the functional form of the bellows 
equations differs only in the algebraic sign of the excitation term, 
which reflects the fact that one bellows is compressed while the 
other is expanded. The response of the turbopump inlet pressure, 

P 6 , in the presence of the velocity excitation, V 4 , is to be deter- 
mined. The component equations are combined manually to produce 
the system equation: 


’ Ps 

Pe/Zj 


= E%DiD3D2Di 


'Zt Q% 
- Qx . 


V^Db Q - DfeDi C 3 - D b D 4 D 3 C 2 j 


(H.5) 


Observe that this system equation is more complex than the system 
equations in the previous examples. The following information is 
immediately available from the system equation: 


(1) NELM = 5 

(2) JBNUM = 3; therefore JTERM (1) = 2 

JTERM (2) = 3 
JTERM (3) =4 

and 

K (1,1) = 5 
K (1,2) = 4 
K ( 2 , 1 ) = 5 
K (2,2) = 4 
K (2,3) = 3 
K (3,1) = 5 
K (3,2) = 4 
K (3,3) = 3 
K (3,4) = 2 
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tation No. Component No. Component Matrix Equations 


1 N 

OTO 


1 to l*> 

cl cr 


1 <r 
□u O 


a? O' 


■in ro 
CL 0 


r * 

a. cr 


oz'cy 


iO <*' 

ora 


CSJ . H- 

■+-* a> 

IK’S 
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m n ii 
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its 
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(3) BSIGN(l) = +1.0 
BSIGN (2) = -1.0 
BSIGN (3) = -1.0 

(4) ITYPE ( 1) = 1 
ITYPE (2) = 7 
ITYPE (3) = 5 
ITYPE (4) = 7 
ITYPE (5) = 1 

(5) SIGN = +1.0 

For computational purposes, it is assumed that Zi = 0 and 
Zt = R. In addition, there are no dissolved gases in the propellant 
(NGAS = 0). Wall elasticity is to be reflected in the speed of sound 
calculations (NELAST = 1). The entire input package is listed on 
the following page. At this point, the reader should make the 
correspondence between each numerical input quantity and the 
FORTRAN statement governing that input. The response magni- 
tude, shown in the output listing, was normalized manually, and 
the resultant response is illustrated in Figure H-8. 
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INPUT 


1FINAL REPuKT EXAMPLE. SWRI PROJECT NO. 02-2888 
5 3 0 1 

1 ? 5?1 
5 3 4 

5 4 

5 4 3 

5 4 3 3 


+.inonoo+oi 
+. 100000+01 
+. 100000+01 
+.300000+08 
+.408000-05 
+ .150000+05 


-. 100000+01 
+.500000+00 
+. 451 U 00+05 
+.bb0000-01 
+. 000000+00 
+.400000+01 


-. 100000+01 
+. b00000+05 
+. 250000+01 
+.000000+00 
+. 000000+00 


+.500000+00 

+.000000+00 

+.000000+00 

+.000000+00 


+. bOOOOO+OS 
+.347000+05 
+. 199000+08 
+.500000+05 


+ .8000 0 0 + 00 
+. 150000+05 
+.800000+00 
+ . 150000+05 


+. 400000-05 
+.400000+01 
+. 400000 - 0 5 
+.400000+01 


+.175000+01 

+.175000+01 


OUTPUT 


FINAL REPORT EXAMPLE. SWRI PROJECT NO. 03-2888 


OMEGA-RAO/SEC 

FREQ-HZ 

b.3 

1.0 


1.5 

12. b 

2.0 

15,7 

2.5 

18.8 

3.0 

25.0 

3.5 

25.1 

4.0 

28.3 

4.5 

31.4 

5.0 

34 » b 

5.5 

37.7 

5.0 

40,8 

b.5 

44.0 

7.0 

4 7,1 

7.5 

50.3 

8.0 

53.4 

8.5 

55.5 

3.0 

58.7 

8,5 

b5.8 

10.0 

bb.O 

10,5 

bR.l 

11.0 

?2.3 

11.5 

75.4 

15.0 

78,5 

12.5 

81.7 

13.0 

84.8 

13.5 

8R.0 

14.0 

81.1 

14.5 

84.5 

15.0 

87,4 

15.5 

100.5 

lb. n 

103.7 

lb. 5 

10 5,8 

17.0 

1 1 0 . 0 

17.5 


TRANS 

TRANS-DB 

775.852 

133.08 

1515.373 

145.0b 

1730.758 

148.13 

2378,342 

155.48 

35b?.?02 

1 b 1 . 8 4 

4 b 5 2 . 0 8 b 

1 b 8 . 8 0 

7582.501 

L 7 7 . 8 b 

14401.82b 

181.50 

53174.748 

50 1 . 05 

82bS, 208 

1 85.58 

4858.452 

1 b 8 . b 5 

5810,558 

158.82 

1558. b21 

147.04 

b 1 4 , b 1 1 

158.45 

511.52b 

10 7. 0 8 

1022. b51 

138. bO 

1801.53b 

151.01 

5840.050 

158,75 

4578. Ib4 

lb ?.53 

bl72.bb l 

174.5b 

8500.758 

185.54 

15035,358 

185.3b 

31378.853 

507.08 

14178b. 521 

237.54 

38505.082 

511.17 

20482.071 

188.5b 

14383.482 

181.45 

11283. 4bl 

18b, b5 

8457.827 

183.03 

8185,881 

180.50 

7585. b8b 

177,80 

bb5b.541 

175.88 

b 1 0 4 . 7 3 8 

174.34 

5b87.4bb 

175.85 


.288000+03 

. 000000+00 
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Figure H-8. Frequency Response Of The Feedline In Example Problem No. 4 
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Problem No. 5: 


As depicted in Figure H-9, the vertical segment of the pro- 
pellant feedline is undergoing forced changes in line length. The 
structural forces applied at Stations 2 and 3 produce measurable 
velocities which, in general, differ in magnitude and/or phase. 
Normally, these velocities can be related through a transfer func- 
tion, G(s). This function has been assigned a simple complex 
variable form in the computer code: G = Gj,-l-iG 2 . The user may 
wish to alter this form to include a specific type of frequency 
dependence based on experimental observations of velocity excited 
lines. Such a program modification can easily be accomplished. 

The length changes imposed on the vertical line segment do 
not influence the pressure-flow relationships of the horizontal seg- 
ment. 


As in the preceding examples, the first step is to describe 
each component by the appropriate matrix expression. The matrix 
equations can then be combined to yield the overall line equation 
which, for this example, is 


" Ps " 
_P 3 /R_ 


= VzQi 


’Pi" 

-Qi. 


v 2 c 2 


(H.6) 


The effect of the structural velocity at Station 3 is contained in the 
column matrix, C s . In computing the response of P 3 to changes 
in line length, the perturbation pressure at the exit to the fuel tank, 
Pj , has been assumed to be negligible. 


At this point, the reader should have no difficulty in con- 
structing a data input package for this problem. For example, the 
input listing may take the form shown below in which the function, G, 
was taken to be ' 

7 r + 1 \/t ‘ 


The computer input and corresponding output for this package, is 
shown on the following page. The modulus of the transfer function 
relating P 3 and V 2 is shown in Figure H-10 
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INPUT 


1 FORCED CHANGE IN LINE LENGTH SWRI PROJECT OP-2883 
2 10 1 
1 3 


innnoo+oi 

loonoo+oi 

+.500000+00 

+ . b00000 + 02 

+.500000+00 

+. b00000+02 

lonnoo+oi 

+.+bl0P0+05 

+.220000+01 

+.000000+00 

+.3+7000+02 - 

300000+08 

+.bb0000-01 

+.000000+00 

+ . 000 0 0 0 + 00 

+.133000+08 + 

+03000-05 

100000+02 

+.000000+00 
+. +00000+01 

+.oonnon+on 

+ .000 n 00 + 00 

+.500000+02 

3 f 1 0 n 0 0 + 02 

+ . +00000+01 

+.707000+00 

+.707000+00 

+.280000+00 

OUTPUT 

FORCED CHANGE 

IN LINE LENGTH SWRI PROJECT 02-2883 

OMEGa-RAD/SEC 

FREQ-HZ 

TRANS 

TRANS-OH 


b.3 

l.o 

1 7 . 3 b b 

5 7 . n 9 


q.+ 

1.5 

17.83 0 

57.33 


12. b 

2.0 

17.88b 

S 7 # b R 


18.? 

2.5 

18.1+8 

5 7.3? 


1 R . H 

3.0 

18. +2 1 

5 8.2? 


2 2.0 

3.5 

18.712 

bb.bb 


28.1 

+ .0 

19.02+ 

5 8.3 1 



+ .5 

..13. 3 b 2 

5 3 . 2 7 


31 . + 

5 . 0 

13.728 

5 9 # b + 


3 + . b 

5.5 

20,127 

b o * n + 


3?. ? 

h # 0 

2 0 . 5 h 3 

8 0 . + 7 


+ n . 3 

b » 5 

2 1 . 0 + U 

■ b 0 * 9 3 


+ + . 0 

7.0 

2 1 . S. b H 

b 1 # + r> 


+ ?.l 

7.5 

2 2.133 

8 1,34 


5 n . 3 

8.0 

22.7? 1 

8 2.41 


5 3.+ 

8.8 

23. +h? 

8 3.11 


Sh . s 

3.0 

2 + . 2 3 3 

8 3.74 


58.7 

3.5 

25.075 

b+. + + 


b2. 3 

1 n # n 

25.33? 

8 5.1 8 


bb*Q 

1 0 . 5 

27.001 

85.32 


b 9 . 1 

11 . 0 

28.083 

b h . ? 0 


? 2 . 3 

11.5 

23.223 

b 7 . 5 H 


? 8 . + 

1 2 . n 

3 0.+ 0 5 

8 8.23 


7 8.5 

12.5 

31.883 

8 3.03 


8 1 . ? 

1 3.0 

3 2 , 5 7 b 

8 3.87 


8 + .8 

13.5 

3 3.3 3 0 

7 n . 1 3 


8 8.0 

1 + . 0 

33.838 

7 0.3 1 


3 1.1 

1+.5 

33.312 

? n . 1 2 


3 + .2 

15.0 

32.2+2 

8 3.+ ? 


87.+ 

15.5 

30. +35 

88.31 


1 0 0,4 

i8.fi 

28.0+5 

88.88 


10 3,7 

lb. 5 

25.311 

8 + . 8 2 


10 8.8 

1 ? . 0 

2 2 . + 7 3 

8?.25 


l lo.o 

17.5 

18.71+ 

5 3.83 


H-5 1 




. 288000+03 

, 000000+00 











